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ABSTRACT 

A  simple  numerical  model  using  a  modified  Euler's  method  was  developed  to 
model  nonlinear  lattices.  This  model  was  used  to  study  the  properties  of  four 
breather  and  kink  type  solitons  in  the  cutoff  modes  of  a  lattice  of  linearly  coupled 
oscillators  with  a  cubic  nonlinearity.  These  cutoff  mode  solitons  were  shown  to 
correspond  very  well  to  the  theoretical  predictions  of  Larraza  and  Putterman  [1984] 
and  the  experimental  work  of  Denardo  [1990].  In  addition,  a  fifth  soliton  was 
discovered  in  the  upper  cutoff  mode,  which  was  not  anticipated  by  the  theory.  A 
preliminary  analytical  attempt  to  describe  this  soliton  and  to  describe  solitons  in  the 
intermediate  modes,  due  to  Larraza,  Putterman,  and  the  author,  is  presented. 
Additional  numerical  work  on  intermediate  mode  solitons  and  domain  walls  was 
performed.  These  studies  showed  that  kink  solitons  are  ubiquitous,  and  that  they 
appear  to  be  intimately  linked  to  domain  wall  structures.  In  order  to  demonstrate 
the  flexibility  of  the  computer  program  developed,  the  model  was  extended  to 
include  two  dimensional  lattices  and  one  dimensional  lattices  with  nonuniform 
characteristics.  Two  dimensional  breather  and  kink  solitons  are  described.  Finally, 
a  Toda  lattice  was  modeled  and  some  preliminary  results  obtained  in  preparation  for 
future  work. 
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I.  INTRODUCTION 

Recent  years  have  seen  a  great  shift  in  emphasis  in  physics  research  as  the  total 
dominance  of  quantum  mechanical  studies  has  yielded  somewhat  to  the  rapidly 
growing  group  of  disciplines  commonly  grouped  under  the  "nonlinear  physics"  title. 
What  is  meant  is  more  specific  --  the  study  of  classical,  nonlinear,  dynamical  systems. 
Within  this  group  of  disciplines  is  the  field  of  soliton  research,  which  boomed  during 
the  last  three  decades.  This  work  has  yielded  practical  results  in  several  areas,  most 
notably  in  the  development  of  fifth  generation  fiber  optic  communications  systems 
which  will  transmit  soliton  pulses,  greatly  increasing  the  permissible  time-bandwidth 
product  of  a  given  system  (Hasegawa  and  Tappert  [1973]).  Solitons,  which  are 
spatially  localized  nonlinear  wave  packets  so  named  because  they  have  many 
particle-like  properties,  have  became  a  key  part  of  the  current  state  of  the  art  in 
cosmology,  particle  physics,  condensed  matter  physics,  and  hydrodynamics,  to  name 
but  a  few. 

Most  of  the  soliton  work  performed  to  date  has  focused  on  continuous  systems, 
such  as  fluids.  Of  the  little  work  done  on  discrete  solitons,  much  has  been 
concentrated  in  the  study  of  cellular  automata.  Recently,  though,  a  group  of 
researchers  at  UCLA  has  been  conducting  research  into  the  characteristics  of  solitons 
in  discrete  lattices.  Experimental  work  has  been  done  by  Denardo  [1990]  on  a 
simple  realization  of  a  nonlinear  lattice  ~  a  one  dimensional  lattice  of  linearly 


coupled  pendula.  He  verified  the  existence  of  two  of  four  predicted  soliton  types 
in  the  cutoff  modes  of  the  lattice  (the  cutoff  modes  are  those  modes  where  all  of  the 
elements  are  either  in  or  out  of  phase  with  both  of  their  nearest  neighbors). 
Additionally,  he  studied  the  behavior  of  these  solitons  as  he  varied  the  frequency 
and  amplitude  of  the  parametric  drive  system  he  used  to  overcome  the  effects  of 
damping. 

The  purpose  of  the  research  presented  here  was  initially  to  provide  numerical 
verification  of  Denardo's  experimental  results,  in  order  to  check  the  degree  to  which 
his  results  actually  corresponded  to  theoretical  predictions.  It  should  be  noted  that 
his  results  were  largely  qualitative,  since  the  coupling  and  damping  parameters  of  the 
lattice  he  used  could  not  be  measured  directly.  A  simple  computer  model  of  the 
pendulum  lattice  was  developed,  tested,  and  then  used  to  quickly  and  accurately 
verify  the  results  of  Denardo.  In  fact,  the  solitons  in  the  pendulum  lattice  were 
found  to  be  in  excellent  agreement  with  the  theoretical  predictions,  which  were 
based  on  a  weakly  nonlinear  approximation.  This  agreement  persisted  even  outside 
the  limits  of  the  approximations  used  in  developing  the  theory,  suggesting  that  these 
soliton  structures  are  very  robust. 

The  linear  lattice  theory  basics  required  for  this  work  are  presented  in  Chapter 
II,  along  with  a  discussion  of  the  theory  of  parametric  drive  and  its  effect  on  the 
stability  of  lattice  dynamics.  Chapter  HI  presents  the  theoretical  development  upon 
which  Denardo's  work  was  based.  Further,  the  results  of  a  study  of  the  stability  of 
various  conditions  of  the  cutoff  modes  is  presented.   This  study  was  performed  in 


order  to  determine  the  effect  that  discreteness  had  on  the  validity  of  the  well  known 
stability  theory  of  Benjamin  and  Feir  [1966],  which  was  developed  for  continuous 
systems.  Finally,  the  numerical  results  which  validated  Denardo's  experimental  work 
are  presented.  The  existence  of  the  third  (which  was  seen  experimentally  by  Larraza 
et  al.  [1990]  in  a  water  channel  experiment),  and  fourth  predicted  solitons,  and  of 
one  previously  unsuspected  soliton  were  also  shown. 

Having  quickly  achieved  the  initial  goal  of  validating  Denardo's  work,  we  went 
on  to  study  modes  intermediate  in  frequency  and  wavelength  between  the  two  cutoff 
modes.  These  results,  along  with  a  brief  sketch  of  a  theoretical  description  which 
is  still  being  developed,  are  presented  in  Chapter  IV.  Many  different  types  of 
soliton-like  structures  were  identified  and  their  properties  studied  in  this  phase  of 
our  work,  leading  us  to  conclude  that,  in  particular,  kink  solitons  (described  in 
Chapter  HI)  appeared  to  be  ubiquitous,  existing  in  every  mode  we  checked.  Moving 
beyond  simple  soliton-like  structures,  the  existence  of  domain  walls  was  shown  for 
many  different  combinations  of  modes  (these  domain  walls  are  essentially  localized 
boundaries  between  two  different  stable  modes). 

Finally,  in  order  to  extend  the  results  obtained  in  the  various  lattice  modes, 
work  is  presented  in  Chapter  V  on  the  effects  of  nonuniformities  in  the  lattice  on 
the  stability  and  motion  of  the  many  different  solitons  and  soliton-like  structures. 
Preliminary  results  from  the  study  of  a  two-dimensional  model  of  the  same  type  of 
lattice,  and  from  the  study  of  a  one  dimensional  model  of  the  Toda  lattice  are  also 
presented  in  Chapter  V.  The  Toda  lattice  is  a  lattice  in  which  the  coupling  between 


nearest  neighbors  is  exponential,  instead  of  linear.  It  has  been  much  studied  in  the 
two  decades  since  it  was  first  formulated  because  it  has  many  attractive 
mathematical  properties  (not  the  least  of  which  is  that  it  is  completely  integrable). 

By  the  time  the  work  presented  here  was  complete,  the  initial  objective  had 
been  broadened  to  include  not  only  a  study  of  nonlinear  dynamics,  but  also  a 
demonstration  of  the  value  of  highly  interactive  computer  modeling  of  real  physical 
systems.  The  numerical  method  used,  which  was  quite  simple,  is  presented  in 
Appendix  A;  a  short  manual  on  using  the  program  and  the  actual  code  are  presented 
in  Appendices  B  and  C.  This  program,  in  its  final  form,  was  designed  to  provide  a 
fast,  accurate  model  of  physical  systems  which  are  composed  of  oscillating  elements 
(such  as  lattices),  which  was  highly  interactive  and  easily  modified.  The  fact  that  the 
program  is  highly  interactive  makes  it  possible  for  the  researcher  to  go  far  beyond 
the  traditional  method  of  numerical  analysis  wherein  a  set  of  parameters  is  put  into 
the  model  and  the  output  is  then  examined.  An  analogy  can  be  drawn  between  the 
traditional  numerical  analysis  methods,  which  act  essentially  as  filters  that  alter  the 
input  in  a  definite  way.  The  model  used  in  this  thesis  is  analogous  to  the  newer 
adaptive  filters,  which  respond  to  the  dynamics  of  the  system  in  a  way  that  increases 
the  utility  of  the  effort.  The  adaptive  features  are  provided  by  the  user  via  keyboard 
commands  which  are  entered  as  the  model  is  running  and  which  take  immediate 
effect. 

The  contribution  which  it  is  hoped  this  thesis  will  make  is,  therefore,  twofold. 
First,  many  interesting  new  results  were  obtained  pertaining  to  solitons  in  nonlinear 


lattices.  These  results  have  served  as  a  launching  point  for  several  interesting 
theoretical  developments  achieved  in  our  group  and  at  UCLA,  and  as  a  means  of 
both  verifying  and  directing  the  work  of  the  experimentalists  (whether  verification 
or  direction  occurred  depended  on  who  got  there  first).  Second,  a  simple  computer 
model  was  developed  which  is  highly  interactive,  somewhat  novel  in  the  level  of 
freedom  it  gives  the  researcher,  and  easily  modified  to  model  other  physical  systems. 
This  program  is  already  being  used  by  several  additional  students,  each  of  whom  is 
making  his  or  her  own  minor  modifications.  It  is  hoped  that  the  program  will  in  the 
end  have  served  as  a  nucleation  site  for  continuing  innovation  in  interactive  physics 
modeling. 


II.  PRELIMINARY:  THE  LINEAR  LATTICE 

A.      THE  SIMPLE  MASS-SPRING  LATTICE  PROBLEM. 

As  a  starting  point  in  the  study  of  solitons  in  nonlinear  lattices,  we  consider 
first  the  behavior  of  linearized  lattices  in  this  chapter,  starting  with  a  simple  lattice 
of  point  masses  connected  by  massless  springs.  Except  in  Chapter  V,  this  and  all 
other  lattices  will  be  taken  to  consist  of  identical  masses  and  stiffnesses.  Also, 
throughout  this  thesis  periodic  boundary  conditions  are  used,  such  that  the  last 
element  is  coupled  both  to  the  next  to  last  element  and  to  the  first  element;  this  is 
essentially  a  finite  ring  lattice,  except  that  all  effects  of  curvature  of  such  a  ring 
lattice  are  ignored.  If  we  restrict  our  analysis  to  longitudinal  motion,  the  exact 
equation  of  motion  for  the  nth  element  is 

x  -—(x   ,+jc   ,-2x),  ELA.1 

*n  V*n+1       n-1    ***iK» 

tn 

with  s  and  m  being  the  stiffness  of  the  springs  and  mass  of  the  point  masses, 
respectively,  and  with  x  being  the  deviation  from  equilibrium  for  each  element. 
Letting 
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we  get 
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We  assume  an  oscillatory  solution  of  the  form 

x  =Aejikna~'*t),  II.A.4 

where  na  is  the  spatial  variable,  and  a  is  the  spacing  of  lattice  elements.  Substitution 
of  this  solution  into  (II.A.3)  and  solving  for  frequency  yields 
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which  is  the  dispersion  relationship  for  this  linear  lattice.  For  a  finite  and  discrete 
lattice,  k  is  quantized.  Thus,  as  seen  in  Figure  II.  1,  the  dispersion  relation  does  not 
give  a  continuous  curve  such  as  the  one  described  by  II.A.5. 

It  is  often  useful  to  work  with  continuum  limit  approximations  of  discrete 
systems,  so  we  will  derive  here  the  continuum  approximation  of  the  linear  ring 
lattice.  If  we  let  y  =  na  be  the  spatial  variable,  and  we  take  the  limit  of  (II.A.3)  as 


-^ 

,-CO 

Aspersion  Curve  for  n=20  Case 

i 

4 

4 

< 

4 

4 



IT) 

cm 

-CM 
*>    (0 

IT) 

< 

1 

4 

■i 

1 

4 

• 

4 

« 

4 

4 

1 

4 

4 

4 

U    c 

1 
MCDCO^CMt-COCO^CMO 

•                        •                       •                        •                                                 •••• 

»-      i-     t-      t-              oooo 

eBaoio  Aouanbaij  jeaun 

Figure  11.1.  Dispersion  relation  for  a  finite  ring  mass-spring  lattice. 


the  spacing  of  the  lattice  goes  to  zero  (and  the  mass  and  stiffness  go  to  zero  and 
infinity,  respectively),  we  get 

^(Vi^r^A,  ii.a.6 
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Thus  one  can  express  the  equation  of  motion  of  the  linear  ring  lattice  in  the 
continuum  limit  by  a  simple  wave  equation: 

H i 1«q  1I.A.7 
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where 
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Finally,  it  is  worthwhile  working  backwards  for  a  moment  to  obtain  a  result 
that  is  interesting  now  and  useful  later.  If  we  start  by  letting  y  be  a  continuous 
rather  than  a  discrete  variable,  and  consider  the  value  of  x(y)  at  two  points  spaced 
a  small  distance  a  to  the  left  and  right  of  a  particular  point  y0,  we  can  express  them 
as  Taylor  series  expansions: 
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Now,  adding  these  together  and  subtracting  twice  x(y0)  gives  us  an  expression 
identical  to  the  right  hand  side  of  (II.A.1),  with  the  exception  of  the  multiplicative 
constant  s,  if  we  ignore  the  higher  order  even  derivatives.  What  does  this  mean 
physically?  If  we  start  with  a  continuous  system  and  discretize  it  as  indicated,  we  get 
the  same  equation  of  motion  as  our  discrete  lattice,  if  we  assume  only  nearest 
neighbors  interact.  Assuming  infinite  interaction  length,  for  example,  means  that  to 
use  (II.A.1)  to  approximate  a  continuous  system  entails  an  error  of  order  4  in  the 
step  size  a  and  in  derivatives  of  x;  for  harmonic  systems,  the  error  is  of  order  (ka)4, 
which  is  small  when  the  second  derivative  is  finite.  However,  when  the  second 
derivative  vanishes,  these  errors  can  be  significant. 

B.      THE  PENDULUM  LATTICE  AND  ITS  GENERALIZATION. 

One  of  the  simplest  experimental  realizations  of  a  nonlinear  lattice  is  a  lattice 
of  coupled  pendula.  Such  a  lattice  was  used  by  Denardo  [1990]  to  study  solitons  in 
discrete  lattices.  An  idealization  of  the  lattice  he  used  is  shown  in  Figure  II.2.  In 
the  actual  lattice,  coupling  was  accomplished  by  tying  knots  between  the  V  shaped 
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Figure  II.2.  Idealized  pendulum  lattice,  from  Denardo  [1990]. 
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strings  which  supported  adjacent  elements,  and  was  assumed  to  be  approximately 
linear.  For  our  purposes,  we  will  assume  a  completely  linear  torsional  spring 
coupling,  for  ease  of  analysis.  Each  pendulum  bob  is  now  acted  on  by  two  forces  - 
gravity  and  the  coupling  force  from  each  nearest  neighbor.  The  exact  equation  of 
motion  is  given  by 

#6 

-TT—^^i+O-i-^J^^.-a  IIB1 

dt2     mL2 


with 


„  -  \8  II-B.2 


This  equation  can  be  rewritten  with  a  coupling  constant  r  replacing  the  coefficient 
of  the  second  term  and  the  sine  term  expressed  as  an  infinite  series: 

3*0.  2  2 


Wn  _  1        <0, 


-Y(e..,-e.,-2e„Ke„-^V„-^e5„+...,        HJ" 


dt2       v     *     "-1       "'         "     6     "    5! 


which  can  then  be  generalized,  with  higher  order  nonlinearities  dropped,  to 


12 


ert-Y(e/l+1+eli.1-2e„)+o>Je/|.ae3n.  n.B.4 

In  this  equation,  the  sign  of  a  is  left  unspecified;  in  Chapter  III,  it  will  be  seen  to  be 
of  critical  importance,  whereas  the  magnitude  is  not.  Throughout  this  thesis,  we  will 
take  a  to  be  either  ±  1  or  0,  as  desired. 

In  order  to  derive  the  dispersion  relation  for  (II.B.4),  we  set  a  equal  to  zero, 
since  the  dispersion  relation  is  a  linear  concept.  As  before,  we  assume  a  solution  of 
the  form  of  (II.A.4),  and  substitute  into  (II.B.4).   Solving  for  frequency  yields 


O) 


N  o>0+4Ysh^- 


II.B.5 


This  dispersion  relation  is  plotted  in  Figure  II.3.  Again,  the  dispersion  relation  is  a 
discrete  function  of  k,  since  we  are  dealing  with  a  finite  discrete  lattice.  As  we  add 
elements  to  the  lattice,  the  density  of  points  on  the  curve  becomes  greater  and 
greater,  until  the  curve  is  continuous.  This  corresponds  to  an  infinite  lattice. 

Denardo  [1990]  studied  solitons  in  the  two  extreme  modes  represented  in 
Figure  II.3,  referred  to  as  the  upper  and  lower  cutoff  modes.  The  lower  cutoff  mode 
corresponds  to  k  =  0,  or  to  an  infinite  wavelength,  and  is  the  mode  where  all  lattice 
elements  oscillate  exactly  in  phase,  so  that  coupling  is  completely  inoperative.  The 
frequency  of  this  mode  is,  as  indicated  by  II.B.5,  exactly  that  of  a  single  linear 
oscillator.   The  upper  cutoff  mode  corresponds  to  k  =  jr/a,  or  wavelength  equal  to 
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Figure  113.  Dispersion  relation  for  generalized  pendulum  lattice. 
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two  lattice  spacings.  In  this  mode,  adjacent  elements  are  exactly  out  of  phase  with 
each  other.  Coupling  is  maximally  effective  in  this  mode,  and  the  frequency  is  given 
by 

g>?-g>;+4Y.  n.B.6 

Taking  the  continuum  limit  of  the  equation  of  motion  exactly  as  before,  we  get 

^_Y^£+sine-o.  ii.b.7 

dt2       dy2 

Keeping  only  the  leading  order  nonlinearity  and  generalizing,  we  get  the  nonlinear 
Klein-Gordon  equation 


^i-Y^-K^e-ae3.  HB.8 

dt2       dy2 


This  is  the  continuum  equation  which  corresponds  to  the  discrete  equation  modeled 
in  this  thesis.  As  with  the  mass  spring  lattice,  the  model  can  be  viewed  as  a  finite 
element  approximation  of  (II.B.7)  with  only  nearest  neighbor  interactions  allowed, 
or  an  approximation  of  (II.B.6)  which  is  accurate,  where  the  second  spatial  derivative 
does  not  vanish,  to  fourth  order  in  ka. 
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C.      THE  DRIVEN  DAMPED  LATTICE. 

The  experimental  work  which  preceded  this  thesis  was  done  using  a  damped 
driven  pendulum  lattice,  and  indeed  the  use  of  drive  and  dissipation  proved  useful 
in  the  numerical  modeling  of  nonlinear  lattices  as  well.  The  reason  for  the  need  for 
drive  in  the  experimental  lattice  is  that  dissipation  is  unavoidable,  so  if  any  steady 
state  was  to  be  studied,  drive  was  necessary.  In  the  case  of  the  work  presented  here, 
the  primary  utility  of  drive  and  dissipation  arises  from  the  fact  that  energy  radiated 
away  from  solitons  as  they  reach  a  steady  state  profile  is  damped  and  does  not 
return  after  traversing  the  ring  lattice.  When  the  free  case  is  studied,  it  is  often 
difficult  to  recognize  the  presence  or  confirm  the  absence  of  solitons  in  the  presence 
of  high  background  radiated  energy. 

In  the  numerical  model,  drive  and  dissipation  were  incorporated  into  the 
equation  of  motion  by  incorporating  a  linear  damping  term,  with  damping  constant 
6,  and  a  modification  to  the  last  term  by  the  inclusion  of  a  term  representing  drive 
of  amplitude  2rj  and  frequency  2<o: 


-Y^i^-r^^K^^4^^^^-^-       UC1 


dt7 


This  is  the  equation  for  a  system  with  acceleration  drive,  as  opposed  to  a  system  with 
displacement  drive,  where  we  would  have  an  additional  factor  of  g>2  multiplying  the 
cosine  in  the  drive  term.    Equation  (II.C.1),  with  the  coupling  term  and  the  cubic 


16 


nonlinearity  neglected,  is  a  form  of  Mathieu's  Equation,  the  analysis  of  which  is  well 
established  (see,  for  example,  Pippard  [1979],  pp.  289-301,  or  Appendix  I  of  Denardo 
[1990]). 
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m.  THE  NONLINEAR  LATTICE  I:  CUTOFF  MODE  SOLITONS 

A.      BRIEF  OVERVIEW  OF  SOLITON  THEORY. 

Solitons  are  "localized  regions  of  motion  with  finite  energy  that  can  appear  in 
media  that  support  waves"  (Denardo  [1990]).  They  can  be  viewed  either  as  finite 
extent  nonlinear  waves  or  as  "particles"  (hence  their  name,  which  was  coined  by 
Zabusky  and  Kruskal  [1965]).  Solitons  have  several  interesting  properties.  They 
exist  because  of  the  competition  between  dispersive  and  nonlinear  effects. 
Dispersive  effects  tend  to  "spread  out"  the  wave,  shedding  energy  via  linear  radiation 
of  energy  at  various  frequencies.  Nonlinear  effects  tend  to  concentrate  the  energy. 
When  the  medium  in  question  obtains  a  profile  that  matches  a  soliton  solution  to  its 
underlying  differential  equations  of  motion,  these  effects  are  exactly  balanced,  and 
energy  is  trapped  in  a  soliton.  Solitons  can  either  be  propagating  or  stationary; 
whereas  most  of  the  theoretical  and  numerical  work  done  until  recently  has 
concerned  itself  almost  exclusively  with  propagating  solitons,  the  work  here  is 
focused  almost  exclusively  on  the  stationary,  or  standing,  solitons  first  discerned 
experimentally  by  Wu  et  al.  [1984]  and  explained  theoretically  by  Larraza  and 
Putterman  [1984].  Another  important  aspect  of  solitons  is  that  they  collide  elastically 
and  are  frequently  very  long-lived.  These  two  facts,  and  the  very  great  localization 
of  energy  within  a  potential  field  that  solitons  represent,  have  made  them  the  subject 
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of  much  interest  in  fields  as  widely  separated  as  molecular  biology  and  particle 
physics. 

Several  well  established  nonlinear  equations  of  motion  have  been  shown  to 
have  soliton  solutions  (Ablowitz  and  Segur  [1981]).  Of  these,  we  are  concerned  here 
with  only  two.  The  first  of  these  is  the  Nonlinear  Schrodinger  (NLS)  Equation: 


dt       dt,1 


III.A.1 


where,  for  a  surface  wave  of  wavenumber  k  on  a  deep  liquid,  the  parameters  are 


v >0 

2  dk2 


UI.A.2 


v-2u>*2 


III.A.3 


and 


Z-x 1. 

dk 


III.A.4 


A  well-established  soliton  solution  to  the  NLS  equation  is 


A- 


N 


— seek 


\ 


2\ 


III.A.5 


where  the  velocity  v  is  a  free  parameter,  and  where 
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b-± 


(*±+i)  iu.a.6 

[dk    2, 


2Y 

An  important  feature  of  this  solution  is  that  the  amplitude  of  the  envelope  divided 
by  the  characteristic  length  of  the  soliton  is  a  constant,  with  value 


K- 


N 


2y  III.A.7 


This  property  proves  very  useful  in  discriminating  actual  soliton  solutions  and 
solutions  which  look  superficially  like  the  hyperbolic  secant  soliton  but  do  not  obey 
this  rule. 

The  other  standard  nonlinear  evolution  equation  with  soliton  solutions  that  is 
relevant  here  is  the  sine-Gordon  equation, 

^i-c2^+(1)jsin(0).o,  in.A.io 

dt2        dx2 

which  is  derived  from  the  equation  of  motion  of  the  pendulum  lattice  and  in  many 
other  physical  applications  (Dodd  et  al.  [1982]).  A  static  kink  solution  to  equation 
(IH.A.8)  is 
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0(jc)-4tan 


-l 


e  c 


This  solution  constitutes  a  transition  between  two  domains  separated  by  360  degrees, 
similar  to  the  separatrix  motion  of  a  single  pendulum.  Solutions  to  the  nonlinear 
evolution  equations  (III.A.1  and  III.A.8)  have  been  shown  to  strictly  meet  the 
accepted  standards  of  what  a  soliton  is.  That  is,  they  collide  elastically  and  have 
finite  spatial  extent.  As  the  exact  equations  of  motion  of  the  generalized  nonlinear 
lattice  are  explored  in  subsequent  sections,  it  will  become  clear  that  these  are  only 
approximate  solutions  to  the  real  lattice.  No  attempt  has  been  made  in  this  thesis 
to  examine  the  elasticity  of  collisions  between  two  or  more  solitons,  so  the  term  is 
used  somewhat  more  loosely  here  than  in  the  literature.  The  width-amplitude  ratio 
test  has  been  used,  however,  to  limit  the  study  to  solutions  at  least  very  closely 
resembling  solitons.  In  fact,  it  is  a  general  observation  that,  in  any  real  systems, 
which  therefore  have  some  dissipative  effects,  however  small,  none  of  the  exact 
soliton  solutions  that  have  been  discussed  in  the  literature  to  date  actually  remain. 
All  of  those  studies  have  been  in  free  systems,  since  a  dissipative  term  in  the 
nonlinear  evolution  equations  invariably  separates  the  equation  form  the  standard 
forms  which  have  exact  soliton  solutions. 


21 


B.      NLS  THEORY  OF  THE  CUTOFF  MODE  LATTICE. 

We  consider  in  this  thesis  the  generalized  equation  resulting  from  a  linearly 
coupled  lattice  of  nonlinear  oscillators.  The  original  system  in  question  was  the 
model  pendulum  lattice  whose  exact  equation  of  motion  is 


o.-YCe^+e^ej-Mne,, 


2««fl  HI.B.l 


where 


(4-£  HI.B.2 

0    L 


and 


Ta2 
mL2 


III.B.3 


and  T  is  the  torsional  constant  of  the  springs  connecting  adjacent  elements,  a  is  the 
lattice  spacing  (which  we  henceforth  take  to  be  unity,  and  drop  from  the  equation), 
m  is  the  mass  of  the  pendulums  bobs,  and  L  is  the  length  of  each  pendulum.  This 
equation  is  usually  linearized  as 

8.-Y(e..I*e..I-2ej-a>;e..  hum 

However,  in  the  weakly  nonlinear  regime,  the  second  term  in  the  Taylor  expansion 
of  the  sine  function  is  included: 
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e„-Y(e„.1+enl-2e^+^e„-^e3„.  UIB-5 

Generalizing  this  equation  as  was  done  in  Chapter  II,  and  adding  drive  and 
dissipation,  we  get  the  equation  used  in  the  computer  model  used  in  this  thesis  (for 
description  of  numerical  methods  used,  see  Appendix  A): 


S.-r(eJI+1+ell.1-2eJ1)-pell-(Wj+2i|co8(2or))e1I+oe 


2.4.9*irvW>/.^YlA    -nvA3.  III.B.6 


Analytically,  this  equation  is  difficult  to  deal  with.  We  consider  first,  therefore, 
only  the  two  cutoff  modes  (see  Chapter  II),  with  linear  frequencies  g>0  and  »i, 
respectively,  where  it  is  possible  to  make  suitable  approximations  and  arrive  at  an 
analytically  tractable  expression  (in  fact,  the  NLS  equation).  We  consider  the 
uniform  lower  cutoff  state  to  be  modulated  by  an  envelope  that  is  slowly  varying  in 
both  space  and  time.  That  is, 

d(x,t)-A(x,t)e>*t+B(xj)e3i(*t+...+c.c,  m.B.7 

where 

tBWAl  HI.B.8 

and  the  lack  of  even  harmonics  is  due  to  the  cubic  nonlinearity.  Substituting  this 
solution  into  the  cubic  Klein-Gordon  (III.B.6),  dropping  all  higher  harmonics,  and 
neglecting  Att  compared  to  A,  because  A(x,t)  is  slowly  varying,  we  get 
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2ju—-y—+(u20-<o2+juP)A+i)A-3a\A?A,  III.B.9 

dt        dx2 


for  the  lower  cutoff  mode  and 


2/a)-^+Y-^i+(cl)5-a)2+ycl>P)A+TiA-3alAFA,  ffl.B.10 

dt       dx2 


for  the  upper  cutoff  mode,  recalling  that 


(d?-<*;+4Y.  raB11 


Thus  both  the  upper  and  lower  cutoff  modes  can  be  represented  by  the  nonlinear 
Schrodinger  equation,  in  the  continuum  limit,  and  at  the  weakly  nonlinear  level.  The 
next  theoretical  task  is  to  examine  the  stability  of  the  uniform  states  of  these  modes; 
following  that,  the  soliton  solutions  of  the  NLS  equation  will  be  examined  in  detail, 
including  their  stability  characteristics. 

C.      STABILITY  OF  THE  CUTOFF  MODES. 

Following  the  method  of  Stuart  and  Diprima  [1978]  as  extended  by  Denardo 
[1990],  we  use  an  amplitude  equation  to  study  the  stability  of  the  uniform  cutoff 
modes  of  our  generalized  lattice.  Stewart  and  Diprima  showed  that  this  method  is 
equivalent  to  the  somewhat  more  complicated  perturbative  methods  used  by 
Benjamin  and  F  ir  [1967]  and  Eckhaus  [1965].  The  basic  idea  is  to  consider  a  small 
perturbation  in  the  sidebands  of  the  uniform  mode  and  to  determine  whether  the 
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mode  is  stable  when  so  perturbed.  The  essential  result  of  Benjamin  and  Feir,  and 
of  Eckhaus  in  a  slightly  different  but  conceptually  equivalent  problem,  was  that,  for 
certain  wavelength  regimes,  resonances  between  the  sidebands  and  the  mode  caused 
energy  to  be  shed  by  the  mode  to  the  sidebands,  which  grow  exponentially. 

Eckhaus  and  Benjamin  and  Feir  were  concerned  with  fluid  phenomena  when 
they  conducted  their  work.  In  fact,  their  results  were  developed  entirely  in  the  realm 
of  continuum  mechanics.  However,  they  can  be  used  profitably  with  continuum  limit 
approximations  of  discrete  lattices,  such  as  the  Nonlinear  Schrodinger  Equation. 
However,  the  effects  of  finite  lattice  size  on  the  theoretical  thresholds  of  instability 
have  not  been  studied  heretofore.  Later  in  this  section  we  will  consider  the  case  of 
a  two  element  oscillator  lattice  under  the  equations  of  motion  given  above.  It  will 
be  shown  that  one  can  solve  approximately  for  the  threshold  of  instability  in  this 
limiting  case,  and  that  the  result  is  in  fact  different  from  that  predicted  by  the 
continuum  theory  (which  is  not  really  surprising;  rather,  it  would  have  been 
surprising  if  the  results  had  been  identical).  Finally,  numerical  results  for  lattices 
from  two  to  40  elements  in  size  will  give  a  very  clear  relationship  between  the 
continuum  theory  and  the  actual  finite  lattice  results.  There  currently  exists  no 
theoretical  framework  in  which  to  place  these  results,  other  than  to  compare  them 
with  the  two  element  lattice  theory  (with  which  the  agreement  is  remarkable). 

1.      Free  Lattice  Stability. 

We  consider  the  NLS  equation  for  the  lower  cutoff  mode  of  the  free 
lattice: 
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27a,-^-c2^+(^-a)2)A-3alAFA.  IU-C.l 

dt         dx2    V  r 


The  uniform  steady  state  is  A  =  Aq,  a  positive  real  constant,  with  the  relationship 
between  amplitude  and  frequency  being 

aZ-L^I-j).  iu.c.2 

To  investigate  the  stability  of  the  motion,  we  set 

A-i^O+Y),  ni.c.3 

where 

T-T(x,r)  UI.C.4 

and 


Substituting  this  into  the  NLS  equation,  neglecting  higher  order  terms  in  Y,  and  using 
(III.C.3)  to  eliminate  zero  order  terms,  we  get  the  evolution  equation  (Denardo 
[1990]) 
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2/a)Tf-r2Trr-3aAo(T  +  T)-0.  III.C.6 

Imposing  a  small  modulation  of  wavenumber  k,  and  including  both  left  and  right 
traveling  waves,  we  set 


V-ae^^+be*'*^, 


III.C.7 


with  a,  b,  and  Q  complex  constants.      Grouping  terms  from  (III.C.26)  according  to 
the  type  of  exponential  factor,  and  setting  the  coefficient  of  each  to  zero,  we  get 


2a)Q+c2fc2-3aAo  -3aA% 

-3a^n  -2(oQ+c2*2-3a/i; 


o  •  J 


a 
_-0. 


III.C.8 


Setting  the  determinant  of  (III.C.28)  to  zero  and  solving  for  Q  (Denardo  [1990]) 
gives 


Q- 


2o> 


1- 


6aA0 


III.C.9 


Thus,  for  Q  to  be  real  and  the  modulation  to  be  stable,  we  require 
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2^2^  III.C.10 


6aAQ<c  2k  \ 


Thus,  for  the  lower  cutoff  mode,  if  a>0  (softening  lattice),  the  motion  will  be 
unstable  at  a  given  amplitude  to  a  sufficiently  long  modulation  wavelength.  Since 
the  wavelength  is  limited  by  lattice  size  for  finite  lattices,  it  is  clear  that  for  finite 
softening  lattices  there  exist  amplitude  thresholds  below  which  the  lower  cutoff  mode 
is  stable  under  modulation.  For  hardening  lattices,  an  identical  procedure  starting 
from  (III.B.10)  yields 

6aAj>c2*2.  ni.C.ll 

Thus,  for  a  hardening  lattice,  the  free  uniform  state  will  be  unstable  for  a  given 
amplitude  for  modulation  wavelengths  beyond  some  threshold,  exactly  analogously 
to  the  softening  lower  cutoff  mode.  Also,  the  softening  upper  cutoff  mode  is  always 
stable  to  sideband  modulations.  This  result  is  in  fact  the  first  of  many  which  show 
the  marked  symmetry  properties  of  the  cutoff  modes. 

2.      Stability  of  the  Uniform  Damped,  Driven  Cutoff  Lattice. 

Starting  from  (IE.B.9),  and  following  the  method  again  of  Denardo  [1990], 
we  consider  first  the  effects  of  parametric  excitation  on  the  nonlinear  cutoff  modes, 
then  the  stability  of  the  damped  driven  lattice  to  sideband  modulation.  Noting  that 
(m.B.9)  is  equivalent  to  (II.C.4),  with  the  addition  only  of  the  cubic  nonlinearity 
term,  it  is  clear  that  the  analysis  of  excitation  from  rest  is  unchanged  from  that 
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performed  in  Chapter  II,  since,  at  and  just  above  the  state  of  rest,  the  additional 
nonlinear  term  plays  no  role.  However,  as  mentioned  in  that  analysis,  the 
nonlinearity  does,  above  some  amplitude,  begin  to  play  a  role  that  counteracts  the 
exponential  growth  of  oscillation  due  to  the  parametric  excitation.  Therefore,  where 
the  linear  case  grows  without  bound  once  the  excitation  from  rest  threshold  is 
reached,  the  nonlinear  lattice  does  not. 

Consider  the   stability   of  these   damped   driven   lattice   states   when 
modulated  by  sidebands.  We  start  as  before  with 

2j<a—-c2— +(<*20-<»2+j<DV)A+r)<»2A-3a\A?A.  III.C.12 

dt        dx2  ' 

To  investigate  continuum  stability  for  a  modulation  of  wavenumber  k,  we  let 

A-A0e»(UV(x,t)),  III.C.13 


where 


^--L^wvW-P2)-  ID.C.14 


3a 


As  before,  we  retain  only  those  terms  linear  in  T,  which  gives  us  (Denardo  [1990]) 
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dt        dx2 


III.C.15 


with 


h- 


2/ w 


III.C.16 


/--^-fo)J-a)2-ya>p±2v/Ti2-p2CD2 
2ycoL 


III.C.17 


and 


1 
^"2/co 


-P]. 


G)0-G>    +/0)| 


III.C.18 


Letting 


Y-ae^^+te-**1^0, 


III.C.19 


we  get  the  coefficient  matrix  equation 


f-hk2-? 
g 


f-hk2-ii[b\ 


in.c.20 


We  thus  find  that  the  motion  is  stable  if 


and 
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&e(f-hk2)z0, 

]f-hkf*tf, 


III.C.21 


III.C.22 


which,  when  the  values  for  f  and  h  are  substituted  in,  yields 


3aAn — < 
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III.C.23 


This  result,  from  Denardo  [1990],  can  be  expressed  graphically  for  an  acceleration 
drive  system,  as  shown  in  Figure  III.l.  To  use  Figure  III.l,  select  any  point  on  the 
tuning  curve.  Decrease  the  ordinate  by  <?\?/2.  If  the  resultant  point  lies  in  a  cross- 
hatched  region,  then  the  state  is  unstable  to  a  modulation  of  wavenumber  k. 
Similarly,  Figure  III.2  shows  the  stability  region  for  the  upper  cutoff  mode,  which  is 
essentially  the  same  as  the  lower  cutoff  mode,  with  w,  replacing  g)0. 

3.      Effects  of  Finite  Lattice  Size  on  Stability. 

The  results  presented  in  sections  1  and  2  were  based  entirely  on  a 
continuum  theory  based  on  the  work  of  Benjamin  and  Feir  [1967].  However,  in 
experimental  and  numerical  work,  finite  lattice  sizes  are  a  necessity,  even  when 
periodic  boundary  conditions  are  used,  as  they  were  in  all  of  the  work  presented 
here.  Theoretically,  although  it  is  clear  that  finite  lattice  size  leads  to  a  discretization 
of  possible  modulation  wavenumbers,  the  effects  of  finite  lattice  size  on  stability 
remain  incompletely  understood.  In  an  effort  to  remedy  this  situation,  we  consider 
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Fig.  m.l.  Stability  of  the  parametrically  driven  steady  state,  from  Denardo  [1990]. 
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Fig.  III.2.   Stability  of  the  driven  upper  cutoff  mode. 
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the  case  of  a  two  element  lattice,  which  can  be  solved  more  directly  (although  still 
only  approximately)  using  the  equations  of  motion.  This  analysis  will  be  performed 
and  the  results  compared  directly  to  those  of  the  previous  Benjamin-Feir  analysis. 
Then  numerical  work  on  the  subject  will  be  presented,  which  connects  in  a 
continuous  fashion  the  results  of  the  Benjamin-Feir  analysis  and  our  straightforward 
N=2  analysis. 

a.     Two  Oscillator  Lattice  Stability. 

For  the  two  oscillator  lattice,  the  exact  equations  of  motion  are 

01  +  a)Jei-ae|-y(81-e2)  ULCJ2A 

e2+G)Je2-ae2-y(ere2).  m.c.25 

Transforming  to  normal  coordinates  according  to  the  transformation 

5-j(61+82)  III.C.26 


c-4(ere2),  iu.c.27 

2   1 


we  get,  after  substitution  into  the  equations  of  motion, 
and 
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?+a&-|(el+eg-o, 


III.C.28 


c+(oJ+2Y)c-5-(e!-eg-o. 


III.C.29 


Some  manipulation  shows  that 


e,±e2- 


2ttZ2+3Z2) 
2C(C2+3£2) 


III.C30 


This  then  yields,  for  equations  of  motion  in  normal  coordinates, 


i+ull-aV+laC2!, 


ra.c.31 


and 


C+((4+2Y)C-aC3+3a$2t. 


III.C.32 


An  examination  of  the  form  of  these  equation  shows  that  each  of  the 
modes,  which  can  be  considered  as  oscillators,  is  parametrically  driving  the  other. 
We  need  to  determine  the  amplitude  threshold  for  excitation  next.  Consider,  for 
weakly  nonlinear  motion,  the  initial  state 


C-Acosut+{higher  harmonics), 


III.C.33 
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z-o, 


III.C.34 


where 


0  4 


For  a  single  parametrically  driven  oscillator, 

0+(Q2+2Ticos2G>f)e-O,  III.C.36 

the  amplitude  threshold  is  given  approximately  by 

hHq2-«4  iilc-37 

This  is  only  valid  if  the  two  frequencies  are  very  similar,  and  is  often  expressed  as 

ti>2QIQ-qL  m.C.38 

We  neglect  the  cubic  term  in  (HI.C.31)  because  the  amplitude  is  small  initially.  The 
other  nonlinear  term  is 


36 


3aA  2cos2o>r£  -  -A  2(  1  +cos2  ur)  i . 

2 


III.C.39 


There  is,  in  addition  to  the  driving  force,  a  restoring  or  antirestoring  force  which 
alters  the  frequency  of  the  driven  oscillator.   The  new  frequency  is  given  by 


Q2-u20-?-aA2. 


III.C.40 


The  amplitude  of  the  parametric  drive  is  given  by 


3,  „2 
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III.C.41 


Excitation  occurs  if 


||aU2>LS-|a^2j-U+2Y-|a^2j 


III.C.42 


The  alteration  of  the  frequency  of  the  driven  oscillator  is  very  important  since  it  is 
a  greater  change  than  the  nonlinear  bending  of  the  frequency  of  the  driver  oscillator: 
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III.C.43 


This  relation  is  not  satisfied  if  a  >  0  (softening  case),  so  the  softening  upper  cutoff 
lattice  is  stable  (as  expected  from  previous  considerations).  For  a  <  0,  it  follows 
from  (III.C.42)  that  the  threshold  condition  is 


^laL42>Y. 
4 


III.C.44 


Therefore  we  can  define  the  critical  amplitude  to  be 


3lol 


IU.C.45 


A  similar  analysis  shows  that,  for  the  lower  cutoff  mode,  the  motion  is  unstable  for 
amplitudes  greater  than  A,..  Before  we  can  compare  this  critical  amplitude,  which 
applies  only  for  the  N  =  2  case,  to  the  Benjamin-Feir  critical  amplitude  (obtained 
from  UI.C10), 
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we  need  to  express  it  in  terms  of  the  original  coordinates  (Equation  (III.C.45)  is  for 
the  normal  coordinate  system): 


AcJ1-2 


_4_~  m.C.47 

N  3a2Y' 


For  a  coupling  parameter  of  0.25,  the  values  are  1.28255  and  0.816497,  respectively. 
It  is  clear  that  the  difference  is  significant  between  the  continuum  theory  prediction 
and  the  discrete  theory;  any  attempt  to  explain  theoretically  the  deviation  from 
continuum  theory  in  the  finite  lattice  limit  must  account  for  this  difference. 

b.     Numerical  Study  of  Finite  Lattice  Effects. 

A  simple  lattice  model,  described  in  Appendix  A,  was  used  to  explore 
the  behavior  of  the  stability  threshold  for  finite  lattices  in  order  to  explore  more 
deeply  the  relationship  between  Benjamin-Feir  continuum  theory  and  the  exact 
theory  just  presented.  In  each  case,  a  lattice  of  amplitude  A  was  modulated  by  a 
small  (about  10%  peak  to  peak  of  A)  modulation  of  wavelength  equal  to  lattice 
length.  The  long  term  behavior  of  the  system  was  then  monitored.  If  the  initial 
amplitude  A  was  above  the  actual  stability  threshold,  then  the  clean,  single  spatial 
frequency  modulation  grew  and  became  complex  spectrally.  Otherwise,  the  system 
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remained  stable  for  long  (at  least  100  cycles).  Various  amplitudes  were  used,  with 
increasingly  fine  resolution,  until  a  good  estimate  of  A,,  was  available.  All  values 
shown  in  graphs  are  accurate  to  within  less  than  one  percent. 

First  we  consider  the  N=2  lattice  as  a  special,  limiting  case  of  the 
finite  lattices.  Figure  III.3  shows  the  theoretical  curves  of  Ac  as  a  function  of 
coupling,  along  with  the  results  obtained  for  a  hardening  lattice  with  numerical  data 
obtained  using  the  model.  It  is  clear  from  the  numerical  results  that  the  approximate 
theory  is  correct,  which  was  expected.  Extension  of  the  exact  two  element  theory 
to  larger  lattices  has  not  yet  been  attempted,  but  a  numerical  determination  of  the 
onset  of  instability  as  a  function  of  lattice  size  was  conducted.  The  results  are  given 
in  Figure  ITL.4.  The  continuum  theory  clearly  determines  the  threshold  for  lattices 
larger  than  20  elements;  for  smaller  lattices,  a  mechanism  probably  related  to  the 
one  studied  in  the  two  element  case  causes  instability  to  occur  sooner. 

While  verifying  this  theory  and  showing  the  extent  of  deviation  from 
continuum  theory  that  the  smallest  possible  finite  lattice  exhibits,  an  interesting  and 
unexpected  additional  result  was  obtained.  Before  describing  it,  it  is  important  to 
point  out  that  the  instability  predicted  by  the  exact  iheory  in  the  last  section  was  due 
to  a  different  mechanism  than  the  Benjamin-Feir  instability.  Instead  of  sideband 
modulational  instability,  the  exact  theory  merely  applied  well  established  (and 
previously  discussed)  parametric  excitation  stability  theory  to  the  case  where  one 
element  is  considered  to  be  parametrically  driven  by  the  second  element.  In  fact, 
one  might  have  been  tempted  to  presuppose  that  the  results  obtained  using  the 
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parametric  drive  theory  would  be  invalid,  since,  as  the  driven  element's  amplitude 
grows,  it  begins  to  be  stabilized  by  its  nonlinearity.  Additionally,  at  some  point  one 
would  guess  that,  in  a  phenomenon  similar  in  appearance  to  the  well  understood 
linear  energy  beating  phenomenon,  the  driven  element  becomes  the  driving  element. 
In  fact,  the  instability  which  occurs  at  the  threshold  amplitudes 
calculated  by  the  exact  parametric  theory  is  a  weak  instability.  For  amplitudes  just 
slightly  greater  than  the  threshold,  the  driving  element  periodically  decreases  and 
then  increases  in  amplitude,  while  the  driven  element  does  the  opposite.  For 
amplitudes  in  the  lower  third  of  the  band  between  the  two  theoretical  thresholds,  the 
two  elements  never  switch  roles.  Then,  there  is  a  band  where  the  roles  switch,  but 
in  a  chaotic  fashion.  That  is,  the  switching  takes  place  at  unpredictable  intervals. 
At  slightly  greater  amplitudes,  the  switching  becomes  more  periodic,  but  the  detailed 
motion  of  each  element  becomes  more  and  more  chaotic.  The  boundary  between 
the  domains  where  switching  does  and  doesn't  take  place  proves  to  be  very  difficult 
to  pin  down,  as  in  fact  do  most  of  the  similar  boundaries  encountered  in  this  study. 
This  can  be  seen,  for  example,  in  the  time  series  of  Figures  LQ.5  and  III.6.  Each  of 
these  is  a  series  of  measurements  of  the  amplitude  of  the  first  of  two  oscillators  in 
a  lattice  with  coupling  of  0.15.  In  Figure  5,  the  lattice  started  with  amplitude  of 
0.8225;  in  Figure  6,  the  initial  amplitude  was  0.84.  For  reference,  the  lower 
threshold  amplitude  for  this  coupling  is  0.65,  and  the  upper  threshold  is  0.99.  In 
Figure  5,  the  lattice  is  just  inside  the  boundary,  that  is,  it  is  just  inside  the  no 
switching  domain  basin  of  attraction.    The  lattice  makes  a  sudden  flip  across  the 
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boundary,  and  then  flips  back  again.  This  behavior  occurred  in  a  free  system  after 
the  system  had  been  running  for  several  hundred  oscillator  cycles.  In  Figure  6,  the 
lattice  lies  initially  just  outside  the  boundary,  then  flips  back  and  forth  a  couple  of 
times  before  settling  inside  the  boundary  where  it  stayed  for  several  thousand 
periods  of  the  oscillator.  The  entire  transient  took  dozens  of  oscillator  cycles,  but 
only  fifteen  of  the  long  time  cycles  of  the  drive.  Based  on  Figures  5  and  6,  one  can 
conclude  that  the  boundary  lies  somewhere  around  0.83.  A  more  striking  example 
of  a  state  lying  exactly  on  the  boundary  is  given  in  Figure  III.7.  The  element  is  from 
an  N=2  lattice  with  maximum  (0.25)  coupling  and  initial  amplitude  of  1.05,  which 
is  51%  into  the  inter-threshold  band.  The  extreme  irregularity  and  clear  lack  of  a 
preferred  state  of  this  free  system  are  typical  of  highly  chaotic  systems. 

While  watching  the  motion  in  the  time  domain  gives  a  vivid  example 
of  chaotic  motion,  it  is  in  the  frequency  and  phase  domains  that  the  best  careful  look 
at  the  onset  of  chaos  can  be  had.  In  order  to  fully  characterize  the  evolution  of  the 
system's  behavior  between  the  two  theoretical  threshold  amplitudes,  a  study  of 
behavior  at  various  amplitudes  was  conducted  using  frequency  and  phase  domain 
methods  .  Figures  III.8  and  IH.9  give  a  sampling  of  the  spectra  of  one  of  the 
elements  at  various  amplitudes.  It  is  important  to  note  that  these  are  representative 
samples  only,  for  the  spectra  for  most  amplitudes  are  continually  shifting,  especially 
in  the  higher  amplitudes  where  there  is  frequent  role  shifting  between  driver  and 
driven  elements. 
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Figure  111.5.   Time  series  for  N=2  element  with  coupling  =  0.15,  amp  =  0.8225. 
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Figure  III.6.  Time  series,  N=2  lattice  element  with  coupling  =  0.15,  amp  =  0.84 
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7igure  III.7.   Chaotic  N=2  lattice  element  time  series  with  max  coupling 
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Figure  III  .7  (cont.) 
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Examining  the  spectra  in  Figures  III.8  and  111.9,  one  sees  that  the 
element  is  monotonic  just  below  the  lower  threshold  (the  exact  theory  threshold). 
As  one  moves  to  just  above  the  threshold,  several  small  tonals  appear  very  close  to 
the  base  frequency,  and  the  signal  becomes  somewhat  more  broadband  in  character. 
At  20%  into  the  "inter-threshold  band",  the  side  tonals  have  become  more 
pronounced,  and  have  moved  farther  from  the  home  frequency.  Note  also  that  they 
are  mainly  on  the  high  side  of  the  main  tonal  -  for  amplitudes  in  the  weakly 
unstable  band,  this  is  typical  of  the  driven  elements.  The  driving  element  typically 
has  sidebands  on  the  low  side  of  the  main  tonal.  At  51%  into  the  band,  a  great  deal 
of  broadband  noise  is  seen;  motion  in  this  regime  is  chaotic.  The  number  of 
sidebands  is  also  much  greater.  Notice,  in  the  first  example  (peak  amplitude  652) 
that  most  of  the  tonals  are  on  the  low  side,  and  that  the  spectrum  is  somewhat 
cleaner  than  the  other.  This  is  indeed  a  driving  element,  although  the  roles  switch 
in  this  portion  of  the  band.  The  second  spectrum  (peak  amplitude  982)  is  typical  of 
an  element  that  is  beginning  to  switch.  A  relatively  clean  spectrum  becomes  very 
noisy,  or  chaotic,  and  then  it  switches  from  predominantly  low  to  high  tonals,  or  vice 
versa.  This  corresponds  to  the  kind  of  chaotic  transient  seen  in  Figures  5,  6  and  7 
in  the  time  domain,  and  is  in  fact  for  the  same  conditions  as  those  which  gave  the 
markedly  chaotic  time  series  in  Figure  IH.7. 

The  next  pair  of  spectra  are  at  72%  into  the  inter-threshold  band. 
The  first  is  remarkably  clean,  and  corresponds  to  a  driven  element,  even  though 
there  is  an  imbalance  to  the  low  side  of  the  main  tonal.  In  fact,  as  one  moves  higher 
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Figure  III.8.  Spectra,  N=2  lattice  elements  just  below  and  above  lower  threshold. 
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Figure  III.9.  Spectra  for  N  =  2  lattice  elements  at  51%  and  75%  of  unstable  band. 
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into  the  band,  one  sees  the  distinction  in  sideband  distribution  between  driver  and 
driven  element  break  down,  and  one  has  to  rely  primarily  on  peak  amplitude  to 
classify.  The  second  element,  based  on  this,  is  clearly  a  driving  element,  and  it  still 
exhibits  the  predominant  low  end  sidebands.  It  is  also  beginning  to  switch,  as 
indicated  by  the  greater  amount  of  noise. 

The  trends  seen  in  this  small  sampling  of  spectra  in  the  inter- 
threshold  band  are  valid  throughout  the  band.  As  one  moves  from  just  below  the 
lower  threshold  to  just  below  the  upper  threshold,  one  sees  a  gradual  move  from 
monofrequency  response  to  a  response  with  many  sidebands.  The  number  of 
sidebands  grows  as  amplitude  increases,  but  it  does  not  grow  with  time  at  a  given 
amplitude  (this  is  different  than  what  one  would  expect  with  a  Benjamin-Feir  type 
instability).  There  is  a  clear  distinction  low  in  the  band  between  driver  and  driven 
element,  both  in  amplitude  and  sideband  location.  This  distinction  due  to  sideband 
location  breaks  down  as  one  moves  above  the  midpoint  of  the  band,  due  to  the 
element  spending  an  equal  amount  of  time  in  each  of  the  roles.  As  seen  above  in 
Figures  III.5  and  III.6,  an  element  near  the  boundary  spends  large  segments  of  time 
in  one  role  or  the  other,  then  may  switch  in  a  chaotic  burst  of  activity  to  the  other 
mode.  Regular  switching  of  roles  between  driver  and  driven  does  not  begin  until 
amplitude  is  well  above  the  boundary  between  the  two  domains  of  attraction.  At 
first  irregular,  or  chaotic,  it  becomes  regular  as  one  moves  up  in  the  band,  which 
corresponds  to  what  has  already  been  seen  in  the  time  domain. 
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Finally,  looking  at  the  lattice  between  the  two  amplitude  thresholds 
in  the  phase  domain  yields  additional  insight  into  the  behavior  of  this  complex  set 
of  states.  Figures  III.  10  and  HI.  11  show  a  succession  of  phase  diagrams  for  the  same 
lattice  element  shown  in  Figures  III.8  and  HL9  in  the  frequency  domain,  and  Figure 
III.7  in  the  time  domain.  The  values  of  amplitude  used  are  also  identical  to  those 
of  Figures  II1.8  and  III.9,  allowing  direct  comparison.  The  first  phase  diagram, 
corresponding  to  an  amplitude  just  below  the  lower  threshold,  is  very  clean  ~  this 
is  stable,  monochromatic  motion,  as  seen  in  the  spectrum.  In  the  next  diagram, 
taken  for  amplitude  slightly  above  the  threshold,  the  same  underlying  shape  is  clearly 
evident  and  takes  the  majority  of  the  "hits"  of  the  Poincare  section,  but  there  is  now 
a  significant  amount  of  scatter  both  in  and  out  of  the  original  shape.  The  idea  of 
"weak  instability"  is  here  given  a  vivid  graphic  representation. 

The  phase  diagram  for  51%,  which  corresponds  to  the  highly  chaotic 
condition  at  the  boundary  discussed  above,  the  original  shape  is  still  visible,  but  its 
banded  structure  has  broken  down.  The  region  is  more  amorphous,  although  some 
hint  of  the  banded  structure  can  still  be  discerned.  The  amount  of  scatter  is  greatly 
increased,  and  there  is  a  slight  darkening  in  the  inner  part  of  the  ellipse, 
corresponding  to  the  fact  that,  with  switching  of  roles  taking  place,  if  on  an  irregular 
basis,  does  lead  the  element  to  spend  an  above  average  amount  of  time  at  small 
energy  levels.  In  fact,  the  well  defined  outer  boundary  of  the  ellipse  corresponds  to 
the  maximum  energy  state,  that  is,  the  state  wherein  all  of  the  system's  energy  is 
concentrated  for  a  brief  moment  in  this  particular  element.  Points  outside  this  ellipse 
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Figure  III.  10.   Phase  diagram  for  just  below  and  above  lower  threshold. 
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Figure  III.  11.   Phase  diagrams  for  51%  and  75%  of  band. 
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are  not  possible  in  a  free  system  with  the  initial  conditions  given.  The  fact  that  this 
boundary  is  seen  here  and  not  in  the  previous  phase  diagram  is  due  to  the  failure, 
low  in  the  weakly  unstable  band,  of  the  driving  element  to  drive  the  other  until  it 
is  completely  spent  (at  which  point  the  other  element  has  all  or  nearly  all  of  the 
system  energy).  Switching  of  roles  cannot  take  place  until  this  outer  ellipse  is 
reached  at  least  occasionally;  with  the  ellipse  so  well  defined,  the  existence  of 
switching,  even  chaotic,  is  reasonable.  At  75%  of  band,  the  second  diagram  in 
Figure  III.  11,  the  regular  switching  of  the  system  is  readily  apparent  from  the 
existence  of  the  dark  orb  in  the  center  and  the  dark  ellipse  without  structure  that 
exists  halfway  out  to  the  maximum  energy  ellipse.  Thus  the  use  of  phase  diagrams 
corroborates  and  illuminates  what  has  been  learned  from  the  time  and  frequency 
domain  results. 

Figures  III.  12,  III.  13  and  III.  14  show  the  behavior  of  lattice  elements 
just  above  the  Benjamin-Feir  threshold.  In  Figure  HI.  12,  the  time  series  of  the  first 
element  of  a  lattice  with  0.15  coupling  is  seen  to  exhibit  regular  switching,  but  the 
period  of  the  switching  is  variable  about  an  average  value.  When  compared  to  the 
series  in  Figure  HI.5,  which  used  an  identical  time  step,  it  can  be  seen  that  the 
average  period  is  roughly  the  same;  in  fact,  this  value  appears  to  be  a  function  of 
coupling  only.  One  could  visualize  the  envelope  of  this  regular  switching  behavior 
as  simply  a  dnoidal  wave  (Figure  EQ.5)  with  kinks  between  every  peak,  although  this 
is  rather  a  bold  step  of  imagination.  Figure  III.  13  shows  various  spectra  for  the 
maximum  coupling  case  (the  same  lattice  as  in  Figures  III.3  and  EQ.4).   The  third 
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Figure  111.12.  Time  series,  N  =  2  lattice  just  above  Benjamin-Feir  threshold. 
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spectrum  was  the  last  taken;  the  growth  of  sidebands  is  clearly  evident.  Finally, 
Figure  III.  14  shows  the  very  diffuse  phase  diagram.  One  can  see  readily,  as  in 
Figure  III.  11,  that  switching  is  occurring,  but  there  is  little  else  to  glean  from  the 
phase  diagram. 

This  rather  detailed  look  at  the  band  between  the  approximate  N=2 
theory  threshold  of  instability  and  the  Benjamin-Feir  continuum  theory  of  instability 
support  the  idea  that  the  two  theories  are  complementary.  When  the  first  threshold 
is  reached,  a  weak  instability  begins  to  manifest  itself  due  to  the  passing  of  the 
parametric  excitation  threshold  amplitude.  As  amplitude  increases,  a  chaotic 
boundary  region  is  traversed  wherein  switching  between  driving  and  driven  states 
occurs  irregularly;  this  boundary  region  passes  into  a  region  where  switching  is 
regularly,  although  not  precisely  periodic.  Above  this,  one  reaches  the  Benjamin- 
Feir  threshold,  where  the  mechanism  for  strong  instability  is  the  dumping  of  energy 
into  successive  sidebands.  The  transition  between  the  two  types  of  instability  is  not 
clear  or  dramatic,  but  it  is  apparently  real. 
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Figure  HI.  13.   Spectra  for  element  just  above  Benjamin-Feir  threshold. 
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Fig.  III.  14.   Phase  diagram  for  an  element  just  above  Benjamin-Feir  threshold. 
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D.      THEORY  OF  CUTOFF  MODE  SOLITONS. 

1.      NLS  Solitons  in  the  Cutoff  Lattice. 

Larraza  and  Putterman  [1984]  showed  that,  for  the  static  case,  with 
frequency  below  the  linear  lower  cutoff  frequency  and  a  softening  lattice,  a  breather 
soliton  is  a  solution: 
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When  the  sign  of  the  nonlinearity  is  changed  to  hardening,  a  solution  was  shown  to 
be,  in  the  p=l  limit: 
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(Denardo  [1990]).   When  considering  the  upper  cutoff  case,  Denardo  showed  that 
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is  a  kink  soliton  solution  of  the  NLS  equation  for  a  softening  lattice,  if  g>0  is  replaced 
by  t^  and 
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is  a  breather  for  the  hardening  lattice. 

With  the  exception  of  this  last  solution,  Denardo  [1990]  observed  all  of  the 
soliton  solutions  for  the  NLS  approximation  of  the  discrete  lattice  experimentally. 
However,  for  the  pendulum  model  results,  the  observations  were  only  qualitative, 
since  it  was  not  possible  to  accurately  measure  the  intrinsic  lattice  parameters  of 
damping  and  coupling,  and  the  coupling  itself  was  probably  not  truly  linear.  This  in 
fact  was  the  starting  point  that  motivated  the  work  of  this  thesis  ~  to  verify  the 
experimental  results  numerically  and  compare  them  with  the  theoretical  predictions. 
This  work  is  presented  in  the  next  section. 

Before  we  proceed  to  the  numerical  work,  some  observations  on  the 
symmetry  properties  of  the  theoretical  predictions  and  their  underlying  physical  basis 
need  to  be  made.  The  symmetry  is  dual  for  the  cutoff  modes  ~  there  is  symmetry 
about  zero  in  a,  and  there  is  symmetry  from  upper  to  lower  cutoff.  We  note  by 
examining  the  equations  for  the  various  soliton  solutions  listed  that  the  parameters 
are  identical  for  both  breathers,  with  the  exception  of  substitution  of  frequencies  and 
sign  changes  which  are  necessary  because  of  the  changes  in  a  and  in  position  on  the 
dispersion  curve.  The  same  holds  true  for  the  kink  solutions. 
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E.      NUMERICAL  OBSERVATIONS  OF  SOLITONS  IN  CUTOFF  MODES. 

In  his  experimental  work,  Denardo  [1990]  showed  that  three  of  the  four 
solitons  predicted  by  the  NLS  theory  for  cutoff  lattices  did  indeed  exist.  He  was 
unable  to  confirm  the  existence  of  one,  the  upper  cutoff  hardening  breather.  In 
addition,  he  mapped  the  drive  planes  of  the  solitons  he  observed.  These  mappings 
consisted  of  varying  drive  frequency  and  amplitude  and  determining  the  region  in 
the  plane  so  mapped  in  which  a  given  soliton  was  stable.  His  lattice  work 
necessarily  suffered  from  a  lack  of  quantitative  comparison  to  theory,  however, 
inasmuch  as  the  actual  coupling  and  damping  parameters  of  the  lattice  were  not 
determined.  This  was  the  initial  starting  point  of  the  research  presented  here  —  to 
provide  quantitative  corroboration  of  Denardo's  results  and  to  compare  them  with 
the  NLS  theory.  This  effort  in  turn  led  to  many  additional  areas  of  exploration, 
which  are  detailed  in  later  chapters. 

Numerically,  all  four  of  the  NLS  solitons  were  found  to  exist  in  stable  states. 
An  additional  soliton-like  structure,  a  hardening  upper  cutoff  kink,  was  also 
identified.  In  the  four  NLS  cases,  the  agreement  to  theory  was  found  to  be 
excellent,  with  the  degree  of  departure  from  theory  apparently  diminishing  linearly 
as  the  time  step  is  decreased.  At  a  time  step  of  one  percent  of  a  period,  the  error 
was  less  than  one  percent,  and  it  got  better  from  there  (see  Appendix  A).  This  close 
match  to  theory  also  naturally  led  to  the  width-amplitude  product  being  constant  for 
given  lattice  parameters,  as  indeed  it  was. 
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In  the  detailed  discussion  of  lattice  model  results  in  this  and  later  chapters,  we 
consider  scaled  displacements  such  that  the  magnitude  of  the  nonlinear  coefficient 
is  unity,  with  only  the  sign  varying  from  hardening  to  softening  systems.  The  natural 
frequency  for  a  single  uncoupled  oscillator  is  also  always  set  to  unity,  except  in 
Chapter  V,  where  nonuniformities  are  discussed.  Thus,  the  lattice  in  each  case  is 
characterized  by  coupling  and  damping  constants.  The  drive  is  characterized  by  its 
amplitude  and  frequency.  Thus,  the  parameter  space  is  four  dimensional;  however, 
only  a  small  part  was  investigated  here.  It  is  also  important  to  note  early  in  the 
discussion  of  actual  lattice  results  that  they  are  not  unique,  for  a  given  point  in  four 
dimensional  parameter  space.  For  a  given  point,  there  can  exist  no  solitons,  there 
may  exist  one  or  more  soliton  solutions,  and  the  solitons  that  do  exist  may  exist  in 
groups  or  singly.  Thus  initial  conditions  play  a  critical  role  in  the  ability  or  inability 
to  achieve  a  desired  state.  Also,  since  the  program  used  was  highly  interactive,  it 
was  possible  to  approach  a  point  in  parameter  space  in  several  different  ways,  each 
having  slightly  different  results.  As  an  examination  of  the  tuning  curves  for 
parametrically  driven  oscillators  makes  clear,  for  example,  the  system  is  much  more 
sensitive  in  general  to  drive  frequency  changes  than  to  drive  amplitude  changes. 

To  give  a  detailed  example  of  the  challenge  offered  by  the  complexity  of  the 
problem,  consider  the  attempt  to  reach  a  state  that  is  in  the  upper  left  corner  of  the 
drive  plane  region  of  stability  for  a  structure  ~  that  is,  a  state  where  drive  frequency 
is  low  and  drive  amplitude  high.  If  one  tries  to  reach  this  state  from  a  precursor 
state  that  is  high  in  amplitude  but  at  a  higher  frequency,  the  transition  will  be 
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impractical  or  impossible  (I  do  not  know  which),  since,  for  any  reasonable  frequency 
increment,  the  transient  at  high  drive  amplitude  is  far  too  great,  and  the  system 
invariably  either  blows  up  or  transforms  into  an  unrelated  lattice  structure.  If,  on 
the  other  hand,  one  approached  the  state  by  starting  from  a  higher  frequency,  lower 
amplitude  point  and  moving  frequency  to  the  desired  value  the  transients  will  be 
manageable,  and  one  can  then  easily  increase  the  drive  amplitude  to  achieve  the 
desired  state.  This  simple  example  points  out  the  characteristics  of  this  work  which 
made  it  very  much  like  trying  to  chart  out  an  unexplored  sea.  It  also  emphasized  the 
great  value  of  interactive  numerical  analysis.  Very  few  of  the  results  achieved  in  this 
thesis  would  have  been  possible  if  a  traditional  numerical  routine  where  one  puts  in 
a  parameter  choice  and  gets  an  output,  then  repeats  the  process,  is  used. 

We  begin  our  look  at  cutoff  solitons  with  the  NLS  cases,  shown  in  Figures 
HI.  15  through  III.  18.  The  upper  and  lower  cutoff  solitons  are,  theoretically, 
symmetrical,  as  discussed  before,  so  we  will  analyze  in  detail  only  the  hardening 
cases.  They  offer  a  richer  harvest,  since  they  are  not  constrained  by  a  potential 
energy  barrier,  and  so  can  exist  at  arbitrarily  high  amplitudes.  The  theoretical  plot 
for  each  is  shown,  and  the  good  agreement  is  apparent.  The  evident  difference  at 
the  low  amplitude  elements  in  the  kinks  are  characteristic  of  all  of  the  numerics  - 
for  a  given  time  step,  the  overshoot  error  is  larger  for  smaller  amplitudes.  As  in  the 
larger  amplitudes,  however,  the  smaller  amplitude  elements  can  be  brought 
arbitrarily  close  to  theory  by  shrinking  the  time  increment  sufficiently. 
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Figure  III.  15.  Lower  cutoff  hardening  kink  example. 
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Figure  III.  16.   Lower  cutoff  softening  breather  example. 
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Figure  III.  17.   Softening  upper  cutoff  kink  example. 
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Figure  III.  18.   Hardening  upper  cutoff  kink  example. 
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While  these  solitons  were  somewhat  difficult  to  achieve  initially,  they  are 
remarkably  robust  structures.  In  fact,  in  the  case  of  the  kink,  the  only  way  to  get  rid 
of  it  is  to  annihilate  it  with  an  antikink  (which  was  done  numerically  on  numerous 
occasions),  or  to  destroy  the  cutoff  mode.  When  one  increases  or  decreases  the  drive 
amplitude,  the  kinks  sharpen  or  broaden  accordingly,  while  maintaining  there 
position  and  identity.  However,  there  is  one  way  that  one  can  force  a  significant 
change  in  the  kinks.  As  the  amplitude  is  increased,  the  kink  sharpens  until  it  gets 
to  a  point  where  the  structure  with  a  node  element  is  unstable  to  random 
perturbations.  The  kink  undergoes  a  transition,  and  arrives  at  a  state  where  the 
node  is  between  two  elements,  as  in  Figure  III.  19.  This  is  effectively  a  one  half 
lattice  site  shift  in  the  position  of  the  kink.  This  phenomenon,  which  will  be  seen 
again  in  Chapter  IV,  is  not  well  understood  theoretically;  in  fact,  as  the  studies  in 
Chapter  IV  will  make  clear,  there  is  in  some  circumstances  a  definite  preference  in 
the  discrete  lattice  for  the  masses  to  be  at  certain  definite  points  of  the  underlying 
waveform  (i.e.,  the  nodes  or  the  peaks,  etc.).  An  interesting  phenomenon  of  the 
breather,  noted  experimentally  by  Denardo,  is  the  existence  in  the  drive  plane 
(Figure  111.20)  of  a  region  where  a  quasiperiodic  state  exists  (Denardo  [1990]). 
Further,  if  one  proceeds  through  this  quasiperiodic  region,  a  region  where  a  highly 
self -focused  state  exists  is  reached.  This  state  was  characterized  in  the  experimental 
lattice  by  one  element  oscillating  with  an  amplitude  of  about  50  degrees  while  the 
two  adjacent  moved  slightly  and  the  remainder  of  the  lattice  was  at  rest.  This 
behavior  is  remarkable  for  a  system  that  is  driven  by  a  global  parametric  drive  - 
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Figure  III.  19.   Hardening  lower  cutoff  kink  with  offset  node. 
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attempt  was  made  to  map  the  drive  plane  to  the  extent  done  experimentally  by 
Denardo  (work  which  will  be  done  by  a  follow  on  student),  the  existence  of  these 
two  states  was  verified  quite  early  on  numerically.  Figure  111.21  shows  the  self- 
focused  hardening  breather.  Since  this  lattice  is  numerical  and  not  subject  to  the 
realities  of  physical  systems,  the  ratio  of  amplitudes  between  the  main  and  adjacent 
elements  is  remarkably  greater  than  1000! 

The  quasiperiodic  state  is  interesting  in  its  own  right,  of  course.  The  period  is 
observed  to  be  stable  in  these  systems;  the  amplitude  of  the  large  element  grows  and 
decays  in  a  regular  cycle.  These  observations  correspond  to  experimental 
observations  made  earlier  by  Denardo  [1990].  This  quasiperiodic  state  can  be 
understood,  I  believe,  as  a  transition  zone  similar  in  general  concept  (although  not 
in  detail)  to  the  transition  zone  between  kinks  with  and  without  nodes  at  masses 
(there  is  a  range  of  parameters  where  either  state  can  exist;  it  is  not  until  one  leaves 
this  range  that  the  transition  takes  place  from  one  to  the  other).  To  the  right  and 
below  the  quasiperiodic  region  in  the  drive  plane  ,  the  breather  essentially  occupies 
five  lattice  sites.  Above  and  to  the  left  of  the  quasiperiodic  region,  it  occupies  only 
three  sites.  In  the  quasiperiodic  region,  it  seems  that  the  two  elements  adjacent  to 
the  high  amplitude  element  first  drive  the  fourth  and  fifth  elements,  and  are  in  turn 
driven  by  the  main  element.  At  some  point,  they  cease  to  drive  the  outer  elements, 
which  then  drive  the  inner  elements  until  they  decay  and  the  process  repeats.  Thus 
the  quasiperiodicity  may  be  viewed  as  a  nonlinear  beating  process  similar  in 
character  to  the  case  studied  in  the  two  element  lattice  in  Section  3.   These  self- 
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Figure  111.20.   Drive  plane  for  upper  cutoff  breather  (from  Denardo  [1990]) 
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links,  and  may  in  some  similar  form  be  significant  as  biological  solitons,  although 
these  connections  are  speculative. 

The  fifth  type  of  cutoff  soliton  was  discovered  quite  by  accident.  It  is  a 
"positive  energy"  hardening  upper  cutoff  kink,  shown  in  Figure  111.22.  By  positive 
energy,  we  refer  to  the  fact  that  the  total  energy  of  a  lattice  with  this  kink  is  greater 
than  that  of  a  uniform  lattice  (the  NLS  kinks  are  by  contrast  negative  energy  kinks; 
in  fact,  the  ability  exists  with  the  lattice  program  to  measure  the  energy  of  kinks  and 
other  structures.  This  was  not  done  due  to  time  constraints,  and  may  be  an  item  of 
interest  in  follow  on  work).  Due  to  the  late  date  of  their  discovery,  only  preliminary 
study  of  them  has  been  conducted.  They  are  effectively  "brightons",  which  are 
solutions  of  the  form 


A(x)-K^l  +asech\Ux-x())  ffl.E.1 


where  K  and  Z  are  dependent  on  lattice  parameters  (Larraza  and  Putterman  [1991]). 
Exact  expressions  for  K  and  Z  have  not  yet  been  derived;  this  theory  is  for 
conceptual  purposes.  Additionally,  these  brightons  appear  to  involve  violations  of 
the  Benjamin-Feir  stability  criteria;  this  will  be  explored  in  more  detail  later. 

The  brightons,  like  the  rest  of  the  cutoff  mode  solitons,  are  remarkably  robust 
structures.  They  maintained  their  identity  even  when  the  time  step  was  increased 
to  25  percent  of  a  period  and  a  perturbation  of  five  percent  was  added. 
Topologically,  the  smallest  brighton  is  one  in  which  exactly  two  elements  are 
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positioning  occurred  at  lower  amplitudes.  As  the  drive  amplitude  was  lowered  to 
0.1,  the  structure  of  Figure  111.22  broke  up  and  underwent  a  prolonged  and  confused 
transient,  ending  up  in  the  state  shown  in  Figure  111.23.  This  is  a  very  complex  state. 
The  brighton  appears  to  have  become  a  darkon,  but  is  equally  clearly  not  an  NLS 
hyperbolic  tangent  kink.  The  mode  itself  has  undergone  a  transition  into  something 
which  preserves  wavelength  at  two  but  is  otherwise  quite  different.  It  is  not  merely 
shifted  spatially,  for  no  shift  could  give  the  amplitude  division  shown  and  still  have 
wavelength  two.  The  only  way  to  describe  it  is  to  think  of  it  as  an  upper  cutoff 
mode  that  has  a  DC  offset,  which  may  also  explain  the  transition  of  brighton  to 
darkon.   Why  this  should  be  so  is  not  understood  even  slightly  at  this  time. 

It  remains  to  be  seen  whether  an  analogous  structure  exists  in  the  softening 
lower  cutoff  mode;  in  fact,  there  may  exist  many  more  solitons  in  the  "simple"  cutoff 
modes.  This  complexity  in  the  simplest  of  modes,  which  can  be  described  by  a 
relatively  simple  evolution  equation  in  which  all  but  the  second  spatial  derivatives 
can  be  ignored  should  lead  one  to  expect  a  far  more  difficult  time  dealing  with 
intermediate  modes,  the  subject  of  the  next  chapter. 
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Figure  111.24.   Final  darkon  like  state  after  brighton  decay. 
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IV.  THE  NONLINEAR  LATTICE  II:  INTERMEDIATE  MODES 

A.      THE  "LAMBDA  FOUR"  MODES. 

The  numerical  study  of  solitons  in  intermediate  modes  requires  first  that  the 
modes  themselves  be  understood.  It  turns  out,  at  least  for  the  "lambda  four"  modes 
(those  whose  wavelength  is  four),  this  is  not  a  trivial  matter.  The  plural  is  used 
because  there  exist  a  variety  of  manifestations,  each  of  which  has  its  own 
characteristics  and  preferred  region  in  the  drive  plane.  The  lambda  four  case  was 
chosen  to  follow  the  cutoff  modes  in  order  of  study  because  it  had  been  much 
observed  experimentally.  The  mode  has  a  natural  linear  frequency  given  by 


«4-y5^27,  WAA 


There  are  three  known  stable  states  for  the  lambda  four  pure  mode,  shown  in 
Figures  IV.  1,  IV.2,  and  IV.3.  During  the  study  of  these  modes,  they  each  acquired 
a  descriptive  vernacular  name.  The  first  to  be  studied  was  the  "plus  zero  minus 
zero"  or  "  +  0-0"  mode,  shown  in  Figure  IV.  1.  It  was  during  the  study  of  this  mode 
numerically,  and  in  parallel  on  the  experimental  lattice,  that  the  preferential 
existence  of  the  "plus  plus  minus  minus"  or  "  +  +  — "  mode  in  certain  regions  of  the 
drive  plane  was  discovered.  The  third  stable  version  (there  are  presumably  many 
more...)  was  discovered  during  the  numerical  study  of  the  "  +  0-0"  mode  drive  plane. 
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The  drive  plane  for  the  "  +  0-0"  state  of  the  lambda  four  mode  was  studied 
extensively  as  a  prelude  to  the  study  of  kink  drive  planes  (the  drive  planes  for  the 
other  two  states  have  not  been  studied  systematically  yet).  Figure  IV.4  shows  the 
drive  plane  study  results  for  one  set  of  initial  conditions.  The  lattice  used  was  of 
forty  elements.  It  is  probable  that,  for  a  lattice  of  different  size,  the  exact  positions 
of  the  lines  shown  may  vary,  but  how  or  whether  they  actually  vary  is  unknown. 
The  very  linear  nature  of  several  of  the  boundaries  of  the  drive  plane  region  of 
stability  was  quite  a  surprise;  the  reason  for  the  linearity  is  not  understood. 

The  behavior  of  the  lattice  at  each  of  the  boundaries  of  the  drive  plane  region 
of  stability  proved  to  be  very  interesting,  and  provided  much  evidence  to  aid  in  the 
understanding  of  the  lattice's  characteristics  in  general  that  proved  valuable  when  we 
moved  on  to  the  study  of  kinks. 

The  behavior  of  the  lattice  at  the  lower  boundaries  proved  to  be  the  most 
complex.  For  all  drive  frequencies  less  than  0.98,  the  lattice  went  unstable  when  it 
reached  the  boundary  shown.  The  instability  was  slow  to  develop.  This  instability 
is  caused  by  the  commencement,  when  amplitude  is  lowered  to  the  curve,  of  very 
slow  and  growing  motion  of  the  line  of  nodes  (the  0's  in  "  +  0-0").  This  is  also  the 
change  that  characterizes  lattice  behavior  when  the  lower  limits  of  the  region  of 
stability  of  the  uniform  mode  are  reached  for  drive  frequencies  greater  than  0.98; 
however,  the  instability  does  not  grow  without  bound  there  and  will  be  discussed 
separately. 
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That  this  in  phase  motion  of  nodes  begins  to  occur  and  that  it  grows  can  be 
understood  qualitatively  by  considering  the  potential  energy  in  which  the  elements 
rest.  In  the  region  of  stability,  the  nodes  are  stable  because  they  rest  at  the  bottom 
of  a  potential  energy  well  caused  by  the  coupling  to  the  adjacent  high  amplitude 
elements.  In  the  pure  mode,  these  adjacent  elements  have  identical  amplitudes,  so 
that  the  two  opposite  coupling  forces  exactly  cancel.  If  we  consider  a  slight 
perturbation  of  the  position  of  one  of  these  nodes,  it  will  feel  less  force  from  the 
element  in  whose  direction  it  moved,  since  the  relative  displacement  has  decreased. 
Conversely,  it  feels  a  greater  force  from  the  other  element.  The  net  force  acting  on 
the  node  element  when  it  is  displaced  differentially  from  its  rest  position  is  then 
restoring,  and  the  nodes  remain  stable.  In  fact,  the  "numerical  temperature"  of  the 
lattice  could  be  determined  by  zooming  in  the  amplitude  scale  until  the  nodes  were 
seen  to  move  --  a  zoom  of  about  240  (or  1012,  the  double  precision  limit!)  The 
motion  is  random,  and  corresponds  by  analogy  to  the  thermal  motion  of  crystal 
elements  in  their  lattice  positions.  Now,  as  the  drive  amplitude  is  lowered,  the 
parametric  excitation  of  the  nodes  decreases.  At  the  same  time,  however,  the 
potential  energy  well  in  which  they  rest  becomes  shorter,  since  the  amplitude  of  the 
antinodal  elements  is  decreasing.  Proceeding  by  analogy,  considering  the  natural 
frequency  of  the  node  if  it  were  to  oscillate  in  the  potential  energy  well  as  a  simple 
harmonic  oscillator,  this  change  in  well  height  would  lead  to  a  decrease  in  natural 
frequency.  Using  equation  (I.C.13)  for  the  threshold  condition  of  parametric  drive, 
one  sees  that  a  decrease  in  natural  frequency  leads  to  a  quadratic  lowering  of  the 
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drive  threshold.  While  not  rigorous,  this  discussion  suggests  how  the  lowering  of  the 
well  leads  to  a  lowering  of  the  drive  threshold  at  a  rate  faster  than  that  at  which  the 
drive  is  lowered.  Thus,  at  some  drive  amplitude,  the  threshold  condition  is  met  and 
node  growth  occurs. 

The  reason  the  node  growth  continues  until  it  destroys  the  mode  (since  this  is 
a  softening  system,  this  merely  means  the  lattice  has  at  least  one  site  which  goes 
over  the  potential  energy  hump)  is  not  well  understood,  although  it  may  be  related 
to  the  Benjamin-Feir  instability  discussed  in  Chapter  II  for  cutoff  lattices.  When  the 
nodes  grow  sufficiently  to  become  commensurate  with  the  amplitude  of  the  antinodal 
elements,  they  continue  to  act  initially  as  a  unified  group.  Thinking  of  them  as  a  sort 
of  lower  cutoff  lattice  embedded  in  the  upper  cutoff  antinodal  lattice,  it  can  be  seen 
that,  for  amplitudes  greater  than  some  threshold  which  depends  on  coupling  and 
lattice  size,  the  "nodal  lattice"  will  become  unstable  to  sidebands.  For  higher  drive 
frequencies,  the  response  amplitude  is  lower,  so  there  should  be  a  threshold 
frequency  above  which  the  Benjamin-Feir  instability  threshold  is  not  reached.  This 
appears  to  be  what  happens  at  drive  frequency  of  0.98. 

Above  this  threshold  frequency,  the  nodes  grow  until  a  steady  state  is  reached 
such  as  that  shown  in  Figure  IV.5.  In  a  phenomenon  similar  to  the  weak  instability 
which  occurred  for  the  N=2  lattice  below  the  Benjamin-Feir  threshold,  a  weakly 
unstable  path  to  chaos  exists  at  the  second  threshold  line  for  higher  frequencies  (see 
Figure  IV.4).  This  threshold  occurs  for  decreasing  drive  amplitude,  which  is  sensible, 
since  it  is  associated  with  increasing  parametric  excitation  due  to  a  still  decreasing 
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potential  energy  well  from  the  antinodal  elements.  Recall  that,  for  the  two  element 
lattice,  a  threshold  of  weak  instability  was  crossed  which  corresponded  to  the 
exceeding  of  a  parametric  drive  threshold;  in  this  case,  the  parametric  drive 
threshold  is  for  decreasing  amplitude,  and  corresponds  to  the  point  where  node 
motion  begins.  The  chaotic  state  of  the  two  element  lattice  occurred  as  one  moved 
further  from  the  drive  threshold,  as  one  crossed  an  apparently  fractal  boundary 
between  the  switching  and  non-switching  basins  of  attraction;  here,  the  switching 
and  non-switching  analogy  corresponds  to  the  switching  between  nodal  lattice  being 
driven  by  and  driving  the  antinodal  lattice.  Although  this  connection  between 
switching  and  non-switching  is  conjectural  to  some  degree,  since  it  is  not  absolutely 
clear  exactly  what  is  happening  at  the  chaotic  threshold,  it  is  clear  that  as  we  move 
further  from  the  threshold  of  parametric  instability,  a  chaotic  band  begins  at  a 
definite  threshold. 

Throughout  this  small  chaotic  region,  the  nodal  lattice  continues  to  maintain 
its  identity  and  oscillates  at  a  mean  amplitude  well  below  the  amplitude  of  the 
antinodal  lattice.  Its  motion  is  jerky  and  aperiodic,  however,  which  typifies  chaotic 
behavior.  As  in  the  case  of  the  two  element  lattice's  approach  to  chaos,  the 
development  of  the  nodes  as  they  approach  the  chaotic  threshold  can  be  illustrated 
well  using  frequency  and  phase  domain  visualizations.  In  Figure  IV.6,  the  spectra 
of  a  nodal  element  are  shown  for  drive  amplitudes  0.075  and  0.056.  The  first  is  well 
above  the  threshold  for  node  motion  (which  is  0.056);  the  spectrum  is  of  thermal 
noise  in  the  rough  vicinity  of  the  drive  frequency.  The  second  is  taken  early  in  the 
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Figure  IV .5.   Stable  state  resulting  from  +0-0  mode  after  nodal  growth. 
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growth  of  the  nodes.  The  frequencies  present  are  1.10  and  0.87,  which  do  not  match 
the  drive  frequency  of  1.0434,  but  are  the  frequencies  about  which  the  thermal  noise 
centered  in  the  first  spectrum.  The  fact  that  the  main  frequency  is  higher  than  the 
drive  frequency  is  consistent  with  the  potential  energy  well  idea,  since  this  additional 
potential  energy  will  increase  the  frequency  of  the  nodal  elements.  The  growth  of 
the  steady  state  oscillations  of  the  nodal  lattice  as  one  lowers  drive  amplitude  toward 
the  chaos  threshold  can  be  seen  in  Figure  IV.7,  which  shows  phase  portraits  of 
several  nodal  and  antinodal  elements  in  their  steady  state  at  drive  amplitudes  of 
0.053  and  0.051.  The  trend,  which  continues  as  one  lowers  drive  amplitude,  is  for 
the  inner  ellipses,  which  are  the  nodal  elements,  to  grow  outward,  with  the  angle  of 
inclination  of  the  ellipses  unchanged.  At  the  same  time,  the  outer  ellipses,  which  are 
the  antinodal  elements  (note  they  appear  in  opposite  quadrants,  since  the  antinodes 
are  in  an  "upper  cutoff"  arrangement),  fatten  up,  with  their  inner  edges  broadening 
and  moving  inward.  As  in  the  two  element  lattice  case,  this  change  is  due  to  the 
greater  amount  of  energy  transferred  from  antinode  to  node  during  each  cycle;  there 
is  as  yet  no  switching  of  roles. 

When  the  drive  amplitude  reaches  0.045,  the  ellipses  are  just  about  to  touch, 
as  shown  in  Figure  IV.8;  the  ellipses  stretched  a  little  further  from  this  state  and 
then  chaotic  behavior  began,  as  seen  in  the  lower  half  of  Figure  IV.8.  In  a  manner 
similar  to  that  of  the  two  element  lattice,  the  ellipses  show  large  scatter.  The  spectra 
of  this  chaotic  state  are  shown  in  Figure  IV.9  for  a  node  and  an  antinode.  The 
characteristic  spread  spectrum  of  chaotic  motion  is  readily  apparent.   However,  as 
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Fig  IV.6.   Spectra  for  nodal  element,   (a)  Thermal  noise,   (b)  Onset  of  growth. 
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Fig.  IV.7.   Phase  portraits  of  steady  states  for  drive  amplitude  =  0.053  and  0.051. 
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Figures  IV.9  and  IV.  10  show,  this  chaos  turned  out  to  be  transient  in  nature 
(although  it  remained  for  many  periods  of  oscillation).  Figures  IV.  10  and  IV.  11  are 
the  same  as  Figures  IV.8  and  IV.9,  only  they  were  taken  after  the  transient  chaos 
died  out.  Figure  IV.  12  shows  the  continuation  of  the  effects  seen  before  when 
lowering  drive  amplitude;  there  were  no  transient  chaotic  states.  Finally,  at  drive 
amplitude  of  0.042,  steady  state  chaos  set  in.  This  is  seen  in  Figure  IV.  13.  Now 
there  is  clear  penetration  of  the  outer  ellipses  by  the  point  of  the  inner  ellipses, 
which  I  believe  represents  graphically  the  point  at  which  switching  of  the  nodes  from 
driven  to  driving  state  occurs. 

An  even  more  vivid  graphical  example  of  the  route  to  chaos  is  seen  in  Figure 
IV.14,  which  is  a  similar  succession  of  phase  portraits  at  various  drive  amplitudes, 
with  the  drive  frequency  set  at  1.054  (higher  than  in  the  previous  example,  which  is 
evident  also  from  the  lower  threshold  amplitude  for  chaotic  motion).  In  the  third 
portrait  of  Figure  IV.14,  the  lower  point  of  the  inner  ellipse  has  clearly  pierced  the 
region  of  the  left  outer  ellipse.  Figure  IV.  15  shows  the  same  elements  a  few  cycles 
later,  with  the  idea  of  chaos  very  obviously  portrayed.  One  point  worth  making 
about  the  change  in  shape  of  the  outer  ellipses  is  that,  for  all  frequencies,  when 
piercing  occurs,  the  tails  of  the  outer  ellipses  stretch  out  and  bend  down  toward  the 
coordinate  (x)  axis;  this  is  because,  as  switching  now  occurs,  these  elements  drive  the 
nodal  elements  until  they  are  nearly  at  rest  in  the  equilibrium  position  (i.e.,  at  the 
origin).  This  state  was  steady  state  (that  is,  it  remained  chaotic,  and  did  not  settle 
out  into  an  ordered  system);  no  transient  chaotic  states  were  observed  at  a  drive 
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Fig.  IV.8.   Phase  portrait  for  drive  amplitude  0.045  showing  transient  chaos. 
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Fig.  IV.9.   Spectra  for  node  and  antinode  showing  transient  chaos. 
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frequency  of  1.054. 

When  drive  amplitude  was  lowered  still  further,  the  lattice  made  a  transition 
from  chaotic  states  to  linear  combinations  of  normal  modes  at  about  amplitude 
0.035.  This  corresponds  to  the  regime  wherein  the  effects  of  the  nonlinearity  in  the 
equation  of  motion  becomes  insignificant.  At  somewhat  lower  amplitudes  still,  about 
0.025,  the  lattice  motion  died  out,  as  the  parametric  threshold  for  steady  state 
excitation  was  no  longer  exceeded. 

The  upper  right  line  in  Figure  IV.4  forms  the  boundary  of  the  region  of 
stability  defined  by  a  similar  instability  to  that  encountered  at  low  drive  amplitudes. 
When  this  line  is  reached,  the  nodes  begin  to  move,  but  they  are  exactly  180  degrees 
out  of  phase  with  each  of  their  node  neighbors.  The  mechanism  for  this  onset  of 
node  motion,  while  probably  related  to  that  which  drives  the  in  phase  motion  just 
described,  is  not  understood.  In  all  cases,  the  motion  grows  continually  until  the 
nodes  and  antinodes  are  commensurate  in  amplitude.  Depending  on  position  on  the 
line  of  instability,  various  types  of  final  states  occur.  For  drive  frequencies  greater 
than  about  1.04,  the  state  that  tends  to  arise  is  the  ++--  mode  (Figure  IV.2), 
although,  with  just  a  slight  push  beyond  the  line  of  instability,  one  can  achieve  the 
mode  shown  in  Figure  IV.3,  which  may  be  preferred  for  certain  drive  parameters 
(this  has  not  been  checked).  For  drive  frequencies  less  than  1.04,  the  transient 
caused  when  the  nodes  become  commensurate  with  the  antinodes  invariably  causes 
the  lattice  to  blow  up  due  to  the  usual  potential  energy  problem.  What  is  most 
remarkable  about  the  out  of  phase  node  motion  boundary  is  that  it  is  very  linear  all 
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Fig  IV.10.   Phase  portrait  for  drive  amplitude  0.045  after  transient  chaos. 
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Fig  rV.ll.   Spectra  for  node  and  antinode,  after  transient  chaos. 
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Fig  IV.  12.   Phase  portraits  for  amplitudes  0.044  and  0.043. 
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Fig  IV.  13.   Steady  state  chaos  -  spectrum  and  drive  portrait. 
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Fig  IV.  14.  Phase  portraits  showing  approach  to  chaos  at  omega  =  1.054. 
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Fig  IV.15.   Phase  portrait  of  steady  state  chaos  at  omega  = 


=  1.054. 
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the  way  to  drive  frequencies  equal  to  the  upper  cutoff  frequency,  where  it  bends 
upward  slightly.  Why  it  should  be  so  is  one  of  the  many  aspects  of  the  lambda  four 
modes  that  is  not  understood  currently. 

The  final  boundary  of  the  stable  +  0-0  region  in  the  drive  plane,  the  upper  left 
curve,  is  characterized  by  a  totally  unexpected  phenomenon  -  isochronous  motion, 
with  respect  to  the  drive,  and  DC  excitation.  When  this  boundary  is  reached,  the 
nodes  move  rapidly  to  a  fixed  DC  offset,  which  is  matched  by  a  DC  offset  in  the 
antinodes,  and  stays  there.  If  amplitude  is  further  increased,  the  offset  magnitude 
increases,  but  at  the  second  curve,  the  antinodes  become  unstable,  breaking  up  into 
complex  modulated  states  which  quickly  lead  to  the  potential  energy  problem  and 
lattice  breakup.  The  nodes  and  antinodes,  when  they  move  to  the  isochronous  mode, 
may  do  so  with  an  imposed  modulation,  depending  on  initial  conditions  as  the 
boundary  is  approached,  the  modulation  stays  fixed  once  the  steady  state  offset  is 
achieved.  Examples  of  modulated  and  unmodulated  isochronously  excited  lattices 
are  given  in  Figures  IV.  16  and  IV.  17. 

As  is  the  upper  right  boundary,  this  boundary  is  very  straight  for  a  large  range  of 
values.  In  this  case,  for  frequencies  less  than  about  0.91,  the  curve  bends  upward; 
also  as  was  the  case  for  the  upper  right  boundary,  the  reason  for  this  profile  is  not 
understood. 

As  this  detailed  look  at  the  drive  plane  region  of  stability  of  the  plain  +0-0 
mode  makes  plain,  the  lambda  four  modes  themselves,  without  solitons  or  other 
complications,  exhibit  quite  complex  and  interesting  behavior  that  is  not  fully 
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Fig  IV.  16.   Unmodulated  isochronously  excited  lattice. 
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Fig  IV.  17.   Modulated  isochronously  excited  lattice. 
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understood.  The  drive  planes  of  the  other  forms  of  the  lambda  four  lattice  have  not 
been  studied,  either  and  may  present  an  additional  range  of  interesting  phenomena. 
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B.      KINKS  IN  THE  LAMBDA  FOUR  MODES. 

As  one  might  expect,  considering  the  complexity  and  number  of  lambda 
four  pure  modes,  there  exists  a  wide  variety  of  kinks  in  these  modes.  Considering 
just  the  +0-0  mode,  we  can  argue  on  purely  geometrical  grounds  for  the  existence 
of  four  types  of  kinks,  corresponding  in  effect  to  phase  shifts  of  90,180,  270  and  360 
degrees.  The  profiles  they  must  exhibit,  shown  in  Figure  IV.  18,  are  based  on  the 
requirement  that  the  entire  lattice  oscillate  at  the  same  frequency.  The  classification 
numbers  given  in  Figure  IV.  18  will  be  used  throughout  the  folio  wing  to  clearly 
identify  the  various  kinks;  the  profiles  are  based  on  a  softening  lattice,  which  was  the 
case  considered  in  the  pure  mode  drive  plane  analysis  described  in  the  previous 
section.  The  main  emphasis  in  this  section  is  indeed  on  softening  +0-0  kinks, 
although  results  are  also  presented  for  hardening  kinks  and  kinks  in  the  other 
lambda  four  modes. 

All  of  the  kink  types  shown  in  Figure  IV.  18  have  been  shown  numerically  to 
exist.  Representative  examples  are  shown  in  Figures  IV.  19  through  IV.22. 

In  addition,  Figure  IV.23  shows  a  Type  IA  kink,  which  is  an  antisymmetric  version 
of  the  Type  I  shown  in  Figure  IV.  18.  Figure  IV.24  demonstrates  clearly  the 
distinction  between  symmetric  and  antisymmetric  kinks.  Presumably,  antisymmetric 
versions  exist  of  the  other  kinktypes  as  well,  although  this  has  not  been  verified. 
They  are  seen  to  be  in  good  agreement  with  the  expected  general  profiles. 
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Fig.  IV.19.   Numerical  results  for  a  Type  I  Softening  +0-0  Kink. 
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Fig.  IV.20.   Numerical  results  for  a  Type  II  Softening  +0-0  Kink. 
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Fig.  IV.22.   Numerical  results  for  Type  IV  Softening  +0-0  Kink. 
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Fig.  IV.23.   Numerical  results  for  Type  IA  Softening  +0-0  Kink 
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Fig.  IV.24.   Symmetric  and  antisymmetric  +0-0  kinks. 
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In  order  to  make  a  clear  comparison  between  kink  behavior  and  normal  mode 
behavior,  a  drive  plane  study  was  conducted  of  a  Type  I  kink.  Figure  IV.25  shows 
the  results  of  this  study,  with  the  outline  of  the  +  0-0  pure  mode  drive  plane  overlaid 
for  comparison.  Several  striking  results  were  obtained  in  the  course  of  this  drive 
plane  exploration,  as  expected,  given  the  complexity  of  the  modal  drive  plane 
boundaries.  In  particular,  the  presence  of  the  kink  produced  dramatic,  localized 
effects  which  appeared  in  many  cases  to  compete  with  the  effects  in  the  wings,  which 
generally  were  identical  to  the  ones  discussed  in  the  pure  mode  drive  plane  results 
of  the  previous  section.  This  is  sensible,  since  in  the  wings,  that  is,  far  from  the  kink, 
the  lattice  should  be  expected  to  behave  as  a  pure  mode. 

Looking  first  at  the  lower  boundary  of  the  stable  region  for  the  Type  I  kink, 
we  see  that  the  general  behavior  is  similar  for  drive  frequencies  above  0.87.  In 
particular,  however,  we  note  that  the  frequency  at  which  node  growth  stabilizes  has 
increased  from  0.98  to  1.02.  This  is  apparently  due  to  the  action  of  the  kink  in 
driving  the  node  growth,  which  occurs  most  rapidly  in  the  vicinity  of  the  kink. 
Chaotic  motion  begins  at  roughly  the  same  spot,  with  drive  frequency  of  1.06  and 
amplitude  of  0.041.  A  particularly  striking  phenomenon  that  was  observed  was  the 
migration  of  the  kink  to  the  right  side  of  the  lattice  (from  the  middle)  when  the 
drive  amplitude  was  lowered  to  0.037,  and  the  chaotic  motion  was  strong.  Figure 
IV.26  shows  this  kink  after  its  migration. 

Far  to  the  left,  below  drive  frequencies  of  0.87,  however,  the  boundary  bends 
upwards  form  that  observed  for  the  pure  mode.  The  reason  for  this  was  not  initially 
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Fig  IV.25.   Drive  plane  for  Softening  Type  I  +0-0  Kink  and  Pure  Mode. 
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Fig  IV.26.   Test  kink  after  migration. 
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clear,  until  suddenly  at  drive  frequency  0.795,  the  kink  went  through  a  severe 
transient  in  which  it  made  a  transition  from  symmetric  to  antisymmetric  states,  with 
the  latter  being  much  sharper,  that  is,  occupying  two  vice  four  lattice  sites.  It 
became  clear  then  that  the  reason  for  the  raised  boundary  is  that,  as  we  lower  drive 
frequency  and  increase  response  amplitude  (since  it  is  a  softening  system),  the  kink 
tends  to  sharpen  up.  By  analogy  to  the  NLS  kinks  for  which  theoretical  descriptions 
are  well  established,  this  sharpening  was  expected;  what  was  unexpected  was  the 
sharp  transient  between  kink  geometries.  This  being  the  case,  it  appears  that  the 
symmetric  kink  can  survive  even  at  very  low  drive  frequencies,  but  only  at  elevated 
drive  amplitudes  (which  of  course  precludes  node  motion).  If  the  amplitude  is 
lowered,  the  nodes  begin  to  move  so  as  to  bring  the  kink  to  its  more  stable  (or  lower 
energy)  state,  the  antisymmetric  state.  As  a  qualitative  check  for  this,  it  was  verified 
that  the  antisymmetric  kink  was  stable  much  closer  to  the  original  boundary  of  pure 
mode  instability. 

For  drive  frequencies  above  1.06,  which  is  close  to  the  linear  frequency  of  the 
upper  cutoff  mode,  the  line  corresponding  to  out  of  phase  node  growth  follows 
closely  that  of  the  pure  mode,  including  the  bending  upward  in  the  vicinity  of  drive 
frequency  1.09,  which  is  the  linear  upper  cutoff  frequency.  However,  below  drive 
frequency  1.06,  the  kink's  stability  boundary  deviates  sharply  from  the  mode's.  From 
that  point  all  the  way  left  to  the  isochronous  curve  (which  matched  that  of  the  linear 
mode  very  closely),  the  boundary  amplitude  increased  in  a  stepwise  fashion.  This 
very  unusual  and  inexplicable  result  was  checked  carefully  at  the  points  indicated, 
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although  it  is  possible  that  the  boundary  exhibits  curvature  in  between  the  data 
points  shown.  This  deviation  from  the  pure  mode  case  extends  to  the  type  of 
behavior  as  well  as  its  location.  In  the  small  region  above  drive  amplitude  0.11  and 
between  frequencies  of  0.99  and  1.04,  there  exist  stable  states  with  out  of  phase  node 
motion;  a  steady  state  example  is  given  in  Figure  IV.27.  Observation  of  the 
development  of  such  a  steady  state  makes  clear  what  the  mechanism  of  the 
departure  from  the  pure  mode  is.  The  kink  apparently  drives  the  nodes  closest  to 
it,  so  that,  for  sufficient  drive  amplitude,  the  node  motion  is  propagated  down  the 
line  of  nodes  and  eventually  grows.  The  node  motion  is  thus  a  radiative  event,  and 
the  node  motion  is  not  precisely  out  of  phase  as  it  is  when  the  "out  of  phase 
instability"  boundary  is  reached  (that  is,  the  upper  right  boundary  of  the  pure  mode 
drive  plane).  When  the  drive  amplitude  is  increased  at  drive  frequency  1.04  to  just 
under  the  original  out  of  phase  instability  curve,  the  lattice  does  indeed  transition  to 
the  ++~  mode  in  a  complex  (three  kinks)  state  (Figure  IV.28),  suggesting  that 
boundary  still  exists,  but  is  just  not  reached  before  the  kink  drives  the  lattice  out  of 
its  original  steady  state  and  into  another.  For  drive  frequencies  below  0.99,  the 
effect  of  the  kink,  which  is  now  sharper  and  of  greater  amplitude,  in  driving  the 
lattice  is  too  great  for  a  steady  state  to  be  reached  with  the  nodes  still  small 
compared  to  the  antinodes;  the  system  "blows  up". 

Not  all  of  the  kinks  observed  numerically  were  of  the  "pure"  types  shown  in 
Figures  IV.  18  through  IV.23.  Frequently,  after  some  perturbation  or  change  in  initial 
conditions  the  lattice  would,  after  a  transient  of  varying  length,  emerge  with  one  or 
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Fig  IV.27.   Steady  state  upper  cutoff  node  motion  in  Type  I  kink  drive  plane. 
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more  complex  kinks.  An  example  is  shown  in  Figure  IV.29.  Using  mismatch 
between  left  and  right  domains  as  our  classification  guide,  we  see  that  this  is  a  Type 
II  kink  (that  is,  we  find  that  antinodes  in  the  right  hand  domain  are  where  they 
should  be  if  there  had  been  no  kink,  but  they  are  "up"  when  they  should  be  "down"  - 
-  Type  II  kink).  Denardo  has  suggested  that  these  complex  versions  of  simple  kink 
types  be  called  "excited  kinks",  which  is  an  attractive  nomenclature,  since  it  recalls 
the  particle  like  qualities  of  solitons  and  confers  on  these  states  the  characteristics 
of  a  particle  which  is  in  an  excited  energy  state.  Another  way  of  viewing  these 
phenomena,  which  ties  in  with  the  discussion  later  in  this  chapter  of  domain  walls, 
is  as  kinks  with  inclusions,  which  might  be  termed  the  solid  state  analogy  —  viewing 
these  solitons  as  analogous  to  transitions  in  crystal  lattices.  Both  views  suffer  from 
a  lack  of  mathematical  underpinnings,  but  they  offer  complementary  insights  into 
what  is  happening  in  the  lattice. 

In  Chapter  II,  the  symmetry  properties  of  NLS  solitons  were  noted  and 
verified.  It  is  a  matter  of  great  interest  to  determine  the  point  on  the  dispersion 
curve,  if  there  is  one,  about  which  this  symmetry  is  based.  Unfortunately,  this 
question  has  not  been  resolved.  However,  some  tantalizing  hints  have  been  found. 
In  particular,  it  has  been  found  that,  if  one  takes  a  Type  I  softening  kink  of  the  +  0-0 
mode  and  performs  the  following  manipulations  on  it,  it  transforms  into  a  positive 
energy  hardening  kink: 


a  -»  -a 


IV.B.l 
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IV.B.2 


where  we  have  relied  on  the  fact  that  we  set 


a>0*l.  rV.B.3 


The  hardening  kink  thus  created  has  a  span  that  is  similar  to  that  of  the  original 
softening  kink,  ai  seen  in  Figure  IV.30.  In  the  cutoff  modes,  such  a  transformation 
would  convert  from  breather  to  kink,  but  here  it  only  takes  us  from  a  negative 
energy  to  a  positive  energy  kink  ~  the  wings  here  still  have  finite  amplitude.  Thus 
the  relation  of  this  phenomenon  to  the  symmetry  between  breathers  and  kinks  is  not 
confirmed;  in  fact,  no  attempt  has  been  made  yet  to  determine  whether  lambda  four 
breathers  even  exist. 

A  classification  of  hardening  +0-0  kinks  can  undoubtedly  be  made  along  the 
lines  indicated  above  for  the  softening  case;  here  we  have  focused  almost  exclusively 
on  the  softening  case  in  order  to  take  advantage  of  the  experimental  work  going  on 
with  the  softening  pendulum  lattice.  An  additional  type  of  hardening  +  0-0  kink  has 
been  observed,  and  is  shown  in  Figure  IV.31. 

One  phenomenon  that  was  observed  in  a  hardening  lattice  that  has  not  been 
observed  elsewhere  yet,  and  which  has  tremendous  possibilities,  was  a  temporally 
phase  shifted  kink  (Figure  IV.32).  Whereas  all  of  the  previous  nonlinear  structures 
studied  in  this  work  and  its  predecessors  have  relied  on  spatial  phase  shifts  and 
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Fig  IV.29.   Complex  Type  II  +0-0  kink. 
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Fig  IV.32.   Additional  type  of  hardening  +0-0  Kink. 
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amplitude  modulation  to  maintain  monofrequency  response,  which  is  the  bottom  line 
requirement  for  steady  state  to  exist,  this  kink  relied  additionally  on  an  eight  degree 
temporal  phase  lead  in  one  of  the  elements  to  maintain  monofrequency  response. 
There  is  no  reason  to  assume  that  there  isn't  an  entire  range  of  phenomena  relying 
on  similar  phase  shifts,  but  thus  far  only  the  one  example  has  been  observed  (to  my 
knowledge). 

The  extreme  richness  of  behavior  encountered  in  the  +0-0  mode  led  to  a 
conscious  decision  to  focus  on  it  to  the  virtual  exclusion  of  the  other  lambda  four 
modes.  However,  there  is  every  reason  to  expect  a  similar  degree  of  complexity  to 
govern  the  behavior  of  these  modes  as  well.  The  large  number  of  additional  modes 
intermediate  between  the  cutoff  modes  have  been  studied  almost  not  at  all,  although 
many  kinks  have  been  seen  in  them  as  a  by  product  of  the  work  presented  here. 
Specifically,  kinks  have  been  observed  in  modes  with  wavelengths  of  three,  five,  six, 
seven,  and  eleven.  The  phenomenon  of  kinks  appears  to  be  ubiquitous. 

C.      DOMAIN  WALLS  IN  THE  NONLINEAR  LATTICE. 

The  phenomenon  of  domain  walls  is  well  known  from  the  study  of  crystal 
structure,  electromagnetism,  and  superconductivity.  A  domain  wall  in  a  lattice  is  a 
(uaually  sharp)  boundary  between  two  different  domains.  An  example  would  be  a 
boundary  between  the  upper  cutoff  mode  and  the  lambda  four  mode.  Typically,  the 
two  domains  are  independent  of  each  other  and  of  the  wall  beyond  a  very  short 
interaction  distance;  that  is,  domain  walls  are  local  phenomena.   They  have  been 
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observed  experimentally  in  the  pendulum  lattice,  although  systematic  study  only 
began  recently,  in  parallel  with  this  work.  There  is  very  little  theoretical 
understanding  of  the  phenomenon,  however,  especially  when  the  transition  between 
domains  takes  place  over  a  finite  length  scale. 

While  domain  walls  had  been  a  topic  intended  for  research  at  the  end  of  this 
work,  if  at  all,  events  intervened  when  they  appeared  on  their  own  during  an 
investigation  of  lambda  four  kinks.  When  the  positive  energy  (Type  II)  lambda  four 
kink  was  first  demonstrated,  the  form  obtained  was  not  that  given  in  the  previous 
section.  In  fact,  it  was  an  "excited"  state,  or  a  state  with  inclusions,  depending  on 
how  one  chooses  to  view  the  phenomenon.  Figure  IV.33  shows  the  original  state, 
and  Figure  IV.34  shows  what  appeared  when  the  indicated  elements  were  removed 
(by  dumping  the  data  to  a  file  and  then  editing  the  file).  As  guessed,  the  removal 
of  the  indicated  elements  did  not  adversely  affect  the  stability  or  identity  of  the  kink. 
The  remarkable  flatness  of  the  "kink",  as  it  was  initially  viewed,  suggested  strongly 
that  it  might  instead  be  a  domain  wall.  To  test  this  hypothesis,  the  elements 
constituting  the  flat  region  of  the  structure  were  removed  and  the  lattice  restarted. 
Not  too  surprisingly,  a  stable  positive  energy  kink  of  normal  size  resulted.  However, 
it  was  also  possible  to  extend  the  flat  region  by  placing  identical  elements  (i.e., 
elements  with  same  amplitude  and  velocity)  in  the  middle  of  it,  and  then  removing 
the  original  primary  domain.  This  resulted  in  a  stable  negative  energy  kink  of  the 
upper  cutoff  mode! 
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Fig  IV.33.   Lambda  Four  Type  II  Kink  with  inclusions. 


129 


c 

CO 

.*: 

F 

cz 

o 

Q 
T3 

o 

1 

o 

D 

+ 

O 

c 

O) 

o 

c 

•  MM 

16 

c 

> 

o 

CD 

E 

"D 

CD 

(3 

DC 

X 

CD 

< 

CX) 

d 


CO 

d 


— i— 
d 


— i— 

CM 

d 


c 
o 

CO 

o 

Q_ 

CD 
O 

CO 


CM 

d 


o 


CO 

d 


(qje)  apnji|dujv 


Fig  IV.34.   Same  kink  as  Figure  34,  with  indicated  elements  removed. 


130 


The  simple  elegance  of  this  result  strongly  suggests  that  domain  walls  might 
profitably  be  viewed  as  regions  where  a  negative  energy  kink  in  one  domain  and  a 
positive  energy  kink  in  the  other  exactly  match  some  set  of  boundary  conditions, 
such  that  the  pair  of  kinks  is  stable,  and  the  domains  on  each  side  are  unaffected  by 
the  presence  of  the  other.  It  is  interesting  to  realize  that  the  use  of  periodic 
boundary  conditions  in  my  model  made  this  conclusion  more  evident,  since  one 
always  needs  and  even  number  of  domain  walls  in  order  to  meet  the  periodic 
boundary  conditions.  Thus  the  idea  of  either  domain  as  an  inclusion  in  a  kink  of  the 
other  is  visually  evident,  whereas  if  domain  walls  had  been  studied  as  structures  in 
their  own  right  from  the  beginning,  the  connection  might  not  have  become  apparent. 

Since  that  auspicious  beginning,  the  only  significant  result  obtained  concerning 
domain  walls  (excepting  the  reaction  to  media  nonuniformities,  which  will  be 
discussed  in  Chapter  V),  is  that  they  too  are  ubiquitous.  This  fact  should  not  be 
surprising,  since  kinks,  the  building  blocks  (evidently)  of  domain  walls,  are 
ubiquitous.  In  fact,  the  wide  variety  of  domain  walls  observed  to  date,  including  the 
extreme  case  of  upper  cutoff/lower  cutoff  domain  walls,  suggests  that,  for  any 
positive  energy  kink,  there  will  exist  a  domain  wall  solution  with  any  mode  that 
exhibits  a  negative  energy  kink,  whose  amplitude  would  be  higher  than  the  original 
mode,  and  vice  versa.  So,  for  example,  in  the  hardening  case,  one  would  expect  an 
upper  cutoff  positive  energy  kink  would  have  domain  wall  solutions  with  any  other 
mode  which  the  lattice  can  support  (i.e.,  limited  only  by  lattice  size  and  mode 
stability  constraints). 
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The  further  study  of  intermediate  modes  is  obviously  necessary,  and  domain 
walls  are  one  of  the  most  important  areas.  Their  potential  bearing  on  many 
technologically  important  areas  in  solid  state  physics  and  critical  point  phenomena 
demands  a  close  examination.  Moreover,  the  need  for  a  comprehensive  theoretical 
treatment  that  allows  all  of  the  phenomena  observed  to  date  is  critical. 
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V.  EXTENSION  OF  THE  BASIC  RESULTS. 

A.      THE  TWO  DIMENSIONAL  LATTICE. 

As  a  natural  extension  of  the  work  reported  to  this  point,  a  two  dimensional 
model  using  the  same  physical  ideas  was  developed.  Only  some  very  preliminary 
work  has  been  done,  but  enough  has  been  discovered  to  warrant  future  work  in  this 
area.  While  no  theoretical  work  has  been  done  for  the  two  dimensional  case  yet,  it 
will  be  useful  to  speculate  a  little  on  the  likely  course  the  theory  will  follow,  since 
it  allows  an  interpretation  to  be  made  of  the  results  presented  that  at  least  seems 
reasonable. 

The  model  used  still  featured  nearest  neighbor  interactions  only,  but  now  in 
two  dimensions.  Diagonal  interactions  are  specifically  ignored,  however.  The  exact 
equation  of  motion  is  given  by 

V.A 


Here  m  and  n  represent  the  row  and  column  number,  respectively.  We  can  see 
quickly,  by  analogy  to  the  one  dimensional  case,  that  the  "upper  cutoff"  case,  where 
each  element  is  exactly  180  out  of  phase  with  each  of  its  four  neighbors,  has  a 
linear  frequency  given  by 
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Wij-wJ+Sy. 


V.A.2 


It  is  also  clear  that,  in  the  case  where  the  lattice  is  upper  cutoff  along  one  of  the 
axes  of  the  lattice  and  lower  cutoff  along  the  other,  that 


wo.i"wijO"wo+4Y» 


V.A.3 


which  of  course  is  the  linear  upper  cutoff  frequency  for  the  one  dimensional  lattice. 
This  is  obvious,  since  the  lattice  sees  no  coupling  in  the  lower  cutoff  axis'  direction. 
In  fact,  the  elimination  of  diagonal  interactions  effectively  decouples  the  two 
orthogonal  directions  x  and  y,  so  that  the  linear  dispersion  relation  can  be  written 
by  direct  analogy  to  that  of  the  one  dimensional  lattice: 


(i)2-o)0+4Ysini 


(k) 
—  +4Ysin' 

UJ 


[2 


V.A.4 


Considering  the  cutoff  cases,  of  which  there  are  now  four,  we  see  that,  again 
by  analogy  to  previous  work, 


#6 

dt2 


■+c 


dx2  +  dy2) 


+g>L  e-ae3, 


V.A.5 


where  m  and  n  are  either  0  or  1.  Without  for  the  time  being  proceeding  to  formal 
proofs,  we  note  the  similarity  of  this  equation  with  the  one  dimensional  case,  and 
speculate  that  the  system  may  be  represented  by  a  two  dimensional  NLS  equation. 
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Fig.  V.l.  Examples  of  possible  2D  cutoff  mode  NLS  solitons. 
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Accordingly,  it  is  reasonable  to  suppose  that  one  might  find  solitons  in  the  two 
dimensional  cutoff  lattice  resembling  the  NLS  solitons  observed  in  one  dimension. 
It  may  be  that  the  soliton  is  truly  two  dimensional,  or  it  may  be  coupled  to  a  pure 
mode  in  the  orthogonal  direction.  Figure  V.l  shows  some  of  the  structures  one 
might  hope  to  find  in  the  cutoff  two  dimensional  lattice. 

Not  surprisingly,  the  two  dimensional  lattice  displayed  fascinating  behavior 
from  the  very  beginning.  The  initial  conditions  used  for  most  of  the  earliest  work 
consisted  of  sinusoidal  modulations  in  the  x  and  y  directions,  with  exactly  one 
wavelength  spanning  the  x  and  y  directions.  In  all  of  these  cases,  the  lattice 
displayed  behavior  that  could  be  divided  neatly  into  three  time  scales.  On  the  first 
time  scale,  of  a  hundred  or  so  periods,  consisted  of  a  rapid  disordering  of  the  lattice 
as  the  initial  disturbance  radiated  energy  (now  in  two  dimensions,  so  the  radiation 
is  more  complex).  Over  the  next  several  hundred  periods,  in  every  case  the  lattice 
resolved  itself  into  a  small  number  of  domains  of  the  original  modulated  cutoff 
mode,  with  each  of  the  domains  being  180  out  of  phase  with  its  neighbors.  Figure 
V.2  shows  the  end  of  such  a  period  for  one  set  of  initial  conditions.  On  the  final 
time  scale,  a  mechanism  which  is  not  well  understodd  but  acts  in  a  fashion  similar 
to  surface  tension  (see  below)  did  indeed  reduce  the  lattice  to  a  single  cutoff 
domain.  This  time  scale  was  often  very  long,  since  in  some  cases  the  surface  tension, 
or  net  difference  in  force  sensed  by  each  domain,  was  often  quite  small.  Figure  V.3 
shows  the  same  lattice  as  seen  in  Figure  V.2,  late  in  this  evolution. 
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Fig  V.2.  A  2D  lattice  at  the  end  of  the  second  time  scale. 
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Fig.  V.3.  The  same  lattice,  late  in  the  final  time  scale. 
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The  results  of  these  initial  investigations  made  it  clear  that,  if  kinks  were  to 
exist,  they  would  have  to  exist  in  some  geometry  which  gave  neither  domain  (there 
are  effectively  only  two,  since  they  are  180°  out  of  phase)  a  dominant  position.  A 
few  obvious  geometries  are  shown  in  Figure  V.4;  all  of  them  did  indeed  support 
kinks.  Figure  V.5  shows  a  very  striking  example,  which  resembles  two  NLS  negative 
energy  kinks,  one  in  the  x  and  one  in  the  y  directions.  This  result  was  very 
encouraging,  if  not  unexpected,  for  it  strongly  supports  the  notion  that  the  theoretical 
description  of  the  two  dimensional  lattice  will  be  very  similar  to  that  of  the  one 
dimensional  lattice.  Of  course,  this  rosy  picture  would  immediately  fall  apart  if 
diagonal  interactions  were  allowed  in  the  model,  as  cross  derivatives  would  then 
abound  in  the  equations  of  motion! 

Figure  V.6  shows  the  same  kink  as  seen  in  Figure  V.5,  but  with  a  complex 
structure  that  resides  at  each  of  the  intersections  between  the  x  and  y  kinks.  This 
was  actually  the  first  two  dimensional  kink  seen,  and  was  used  as  the  starting  point 
in  getting  to  Figure  V.5.  The  structures  seen  are  very  stable,  and  they  are  certainly 
not  yet  understood. 

The  final  quick  foray  into  the  two  dimensional  arena  which  was  undertaken 
was  an  effort  to  determine  if  there  existed  stable  breathers.  Initially,  it  was  felt  that 
such  structures  might  not  exist,  due  to  the  existence  of  diffraction  effects.  However, 
they  were  found  to  exist;  in  fact,  it  turns  out  that  if  one  starts  with  a  lattice 
completely  at  rest  and  then  kicks  one  element  (preferably  the  middle  one  for 
visualization  purposes)  with  a  reasonable  amplitude  (i.e.,  enough  that  nonlinearities 
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Fig.  V.4.  Possible  geometries  for  stable  2D  kinks,  where  surface  tension  is  zero. 


140 


Fig.  V.5.  An  actual  2D  lattice  with  symmetric  kinks. 
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Fig.  V.6.  The  same  lattice  with  complex  structures  at  the  kink  intersections. 
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can  play  a  significant  role),and  in  fact  a  great  deal  of  radiation  occurs  initially  which 
leaves  behind  a  very  clean  two  dimensional  breather,  seen  in  Figure  V.7.  It  is 
actually  easier  to  get  a  two  dimensional  breather  than  a  one  dimensional  one, 
because  the  radiation  spreads  out  in  two  dimensions  and  is  much  more  rapidly 
attenuated.  Even  in  the  free  case,  a  very  clean  breather  can  be  obtained,  since  the 
energy  radiated  away  as  the  breather  seeks  its  equilibrium  shape  has  to  spread  itself 
over  the  entire  two  dimensional  span  of  the  lattice,  resulting  in  very  low  amplitude 
radiation  passing  through  the  breather  at  any  given  time  (although  of  course  there 
is  always  some). 

A  theoretical  consideration  which  was  suggested  for  the  first  time  by  the  results 
obtained  with  the  two  dimensional  lattice  is  the  idea  of  "surface  tension"  at  domain 
boundaries.  Suppose  for  the  moment  that  two  similar  domains  are  connected  by  a 
kink  of  some  two  dimensional  shape.  An  example  which  will  motivate  the  discussion 
is  given  in  Figure  V.8.  The  question  arises  whether  one  or  the  other  of  the  domains 
will  dominate  the  other;  that  is,  will  one  domain  gradually  force  all  of  the  elements 
of  the  other  to  shift  modes?  As  it  turns  out,  this  is  indeed  the  case.  In  fact,  for  two 
domains  of  similar  amplitude  (i.e.,  identical  domains  separated  by  kinks),  the  domain 
which  is  concave  will  dominate  the  convex  domain,  regardless  of  the  size  of  the 
domain. 

That  this  is  the  case  is  reasonable,  when  the  issue  is  considered  as  analogous 
to  the  phenomenon  of  surface  tension.  Since  the  concave  domain  exerts  more  force 
per  boundary  element  on  the  convex  domain  than  vice  versa  (because  it  has  more 
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Fig  V.7.  2D  breather  resulting  from  kicking  of  middle  element. 
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elements  interacting  at  the  boundary),  there  is  an  effective  surface  tension,  or  net 
force,  acting  on  the  convex  domain.  Since,  unlike  the  case  of  an  air  filled  bubble 
being  compressed  in  water,  there  is  no  reaction  force  which  builds  up  as  compression 
takes  place,  the  convex  domain  is  gradually  "compressed"  away.  It  should  be  noted 
that  this  phenomenon  depends  on  the  fact  that  coupling  of  an  element  to  its  home 
domain  is  negligible,  so  that  one  can  simply  add  up  the  coupling  forces  to  the 
opposite  domain  to  determine  the  dominant  domain. 

Presumably  a  modified  version  of  this  rule  would  be  operative  when  dissimilar 
domains  are  in  contact;  however,  it  is  not  immediately  clear  what  exact  form  the  rule 
would  take  after  the  differences  in  mean  amplitudes  are  accounted  for.  In  any  case, 
no  numerical  data  for  such  cases  has  so  far  been  collected,  so  it  remains  a  truly  open 
question  (as  does  virtually  everything  having  to  do  with  two  dimensional  lattices 
given  by  (V.A.I)!). 

This  brief  extension  of  the  lattice  model  into  two  dimensions  leaves  much  left 
to  explore.  It  is  clear,  though,  that  there  are  many  new  phenomena  likely  to  reward 
such  efforts.  The  model  also  permits  study  of  surfaces,  including  radiation  and  finite 
structures,  providing  that  sufficient  computer  power  is  available  to  provide  a  small 
enough  discretization  of  the  surface  (that  is,  more  memory  is  needed,  since  the 
model  is  limited  on  the  machine  used  to  a  40x40  lattice,  which  is  not  adequate  to  use 
for  a  surface  model). 
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B.      THE  NONUNIFORM  LATTICE. 

As  a  second  attempt  to  extend  the  work  reported  in  Chapters  III  and  IV,  an 
investigation  was  made  into  the  effects  of  nonuniformities  on  lattice  behavior.  To 
some  extent,  this  work  was  motivated  by  results  obtained  on  the  experimental  lattice, 
since  it  was  assumed  that  that  lattice  was  nonuniform,  and  the  effects  and  indeed 
types  of  nonuniformities  needed  to  be  understoodif  the  experimental  work  is  to  be 
joined  to  the  numerical.  Additionally,  it  was  desired  to  test  a  method  of  removing 
the  effects  of  nonuniformities  which  had  been  tried  by  the  group  at  UCLA  with 
which  we  closely  worked.  This  method  consisted  of  first  measuring  the  amplitudes 
of  each  element  in  the  lattice  in  a  pure  cutoff  mode,  then  dividing  these  results  into 
the  measured  lattice  amplitudes  in  the  presence  of  solitons.  While  this  method  was 
based  on  intuition,  it  proved  to  give  smooth  results  from  results  that  had  been 
difficult  to  interpret.  This  method  of  division  greatly  facilitated  their  study  of 
solitons,  but  it  needed  to  be  verified  before  full  confidence  was  placed  in  it. 

Originally,  variations  in  coupling,  natural  frequency,  and  nonlinear  coefficient 
were  tried.  The  nonlinear  coefficient  was  found  to  affect  results  only  slightly. 
Accordingly,  only  results  from  nonuniformities  in  coupling  and  natural  frequency  are 
discussed  in  this  chapter.  Of  these,  natural  frequency  was  the  most  studied,  although 
only  because  of  time  constraints. 

The  parameters  were  varied  in  three  different  methods  ~  randomly,  with  a 
linear  gradient,  and  with  sinusoidal  gradients.  The  program  allowed  the  user  to 
choose  the  type  of  variation  and  the  amount.  Also,  it  was  during  these  investigations 
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that  a  waterfall  display  was  added  to  the  program  so  that  a  parameter  could  be 
instantaneously  changed  by  one  of  the  three  methods  just  mentioned,  and  the  impact 
viewed  in  a  format  which  made  changes  easier  to  identify. 

Finally,  before  proceeding  to  discuss  actual  results,  it  should  be  noted  that  no 
attempt  has  yet  been  made  to  fit  these  and  the  experimental  results  into  a  theoretical 
framework.  This  work  follows  in  the  time  honored  tradition  of  forging  ahead  with 
experimentation  in  order  to  provide  grist  for  the  theoreticians'  mill,  although  I 
frankly  would  have  preferred  to  have  tended  that  mill  myself,  had  time  allowed! 

The  first  result  was  obtained  early  in  this  thesis  research,  using  an  older  version 
of  the  program  in  which  the  nonlinear  coefficient  was  allowed  to  vary.  It  was  found 
that  the  division  method  developed  by  the  UCLA  group  [Putterman,  unpublished, 
1990]  was  valid.  An  example  is  shown  in  Figures  V.9  through  V.ll.  The  rough 
soliton-iike  structure  in  Figure  V.9  was  divided,  element  by  element,  by  the 
amplitudes  of  the  cutoff  mode  in  Figure  V.10  to  produce  the  smooth  soliton  seen  in 
Figure  V.ll.  The  method  was  found  to  be  valid  for  both  coupling  and  natural 
frequency  variations.  One  additional  item  noted,  which  served  as  the  starting  point 
for  later  work  on  nonuniformities,  was  that,  for  softening  nonlinearities,  the  solitons 
invariably  moved  to  the  local  minimum  in  the  varied  parameter.  Since  reduction  in 
coupling  leads  to  a  reduction  in  frequency,  this  result  means  that,  in  general, 
softening  kinks  move  to  the  local  minimum  frequency  in  the  lattice. 

In  order  to  characterize  the  behavior  of  kinks,  breathers,  and  domain  walls  in 
various  types  of  nonuniform  lattices,  we  returned  to  this  line  of  inquiry  recently. 
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Fig.  V.9.   Soliton  in  a  random  coupling  lattice. 
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Fig.  V.10.  Actual  variation  in  coupling  of  lattice  in  Fig.  1,  with  pure  mode  shown. 
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The  initial  results  have  been  consistent,  but  as  many  questions  have  been  raised  as 
have  been  answered.  The  first  result  in  the  second  round  of  nonuniformity  work 
complemented  the  last  result  of  the  first  round  ~  hardening  kinks  move  to  the  local 
maximum  frequency  in  the  lattice.  In  the  case  of  continuous  gradients  in  either 
natural  frequency  or  coupling,  the  hardening  kinks  moved  steadily  to  the  right 
(increasing  frequency  direction)  until  they  reached  the  end  of  the  lattice,  where  there 
was  a  discontinuous  drop  in  frequency  due  to  the  periodic  boundary  conditions. 
Since  there  were  always  two  hardening  kinks,  this  provided  a  convenient  way  to 
collide  kinks,  since  the  leftmost  kink  inevitably  caught  up  with  the  pinned  right  hand 
kink.  As  expected,  these  encounters  produced  annihilations  of  both  kinks,  leaving 
an  undisturbed  nonuniform  lattice  mode. 

In  order  to  better  study  collisions,  and  to  discern  whether  the  motion  was  a 
diffusion  process  or  a  process  of  motion  in  a  potential  field,  a  sinusoidal  variation 
of  coupling  was  introduced,  with  one  kink  on  either  side  of  one  of  the  peaks  in  the 
variation.  In  particular,  we  wanted  to  find  out  whether,  when  a  kink  arrived  at  a 
maximum  frequency,  it  overshot  the  mark  and  exhibited  oscillatory  behavior.  In 
fact,  this  did  not  happen  at  all,  so  we  were  also  unable  to  observe  two  kinks  passing 
through  one  another  in  this  way  (since  they  both  stopped  at  the  local  maximum). 
Nevertheless,  this  sinusoidal  parameter  variation  feature  was  found  to  be  useful  for 
separating  kinks,  since,  if  one  places  a  minimum  between  two  kinks,  they  will  move 
in  opposite  directions.  Since  the  variations  could  be  turned  off  as  easily  as  on,  one 
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could  separate  two  kinks  using  a  sinusoidal  variation  and  then  turn  off  the  variation 
to  observe  the  cleanly  separated  kinks  in  a  uniform  lattice. 

We  turned  next  to  domain  walls,  anxious  to  see  whether  they  moved 
preferentially  in  one  direction  or  the  other.  In  fact,  no  matter  how  hard  we  pushed 
the  lattice,  we  could  not  get  domain  walls  to  move.  The  UCLA  group  independently 
reported  similar  results.  Since  one  possible  mechanism  for  this  pinning  of  domain 
walls  was  that  negative  energy  and  positive  energy  kinks  might  tend  to  move  in 
opposite  directions  (viewing  domain  walls  in  this  instance  as  paired  kinks),  I 
immediately  set  after  positive  energy  kinks,  to  see  if  they  moved.  Unfortunately, 
they  also  appear  to  be  pinned.  Why  this  should  be  so  is  not  clear,  and  the  fact  that 
it  is  so  can  be  interpreted  as  an  argument  against  the  notion  that  kinks  and  domain 
walls  are  intimately  connected  phenomena.  On  the  other  hand,  it  may  simply  be 
that  the  departure  from  normal  mode  amplitude  by  individual  lattice  elements  is  the 
key  cause  of  soliton  motion  in  a  nonuniform  lattice.  If  this  were  the  case,  and  it 
doesn't  seem  unreasonable  to  suppose  that  it  might  be,  domain  walls  and  positive 
energy  kinks  alike  would  not  move,  or  would  move  so  slowly  that  the  experiments 
performed  to  date  would  not  have  detected  the  motion,  since  the  changes  are  small 
in  these  cases.  Another  way  of  expressing  the  idea  is  to  note  that,  in  the  negative 
energy  kinks  considered  here  (only  a  small  subset  of  those  possible),  the  envelope 
one  would  draw  of  the  lattice  often  changes  sign  across  the  kink,  and  always  makes 
a  major  transition  in  the  kink  region.  In  the  positive  energy  kinks  and  domain  walls, 
and  possibly  "darkons",  the  envelope  only  changes  slightly,  and  retains  the  same  sign 
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throughout.  One  possible  way  of  checking  this  idea  is  to  measure  the  speed  of 
motion  of  negative  energy  kinks  as  a  function  of  amplitude,  since  presumably  if  the 
hypothesis  were  valid,  they  would  move  more  rapidly  as  amplitude  increased.  In  any 
case,  this  area  is  one  which  might  fully  occupy  a  future  thesis  student. 

For  certain  types  of  abrupt  changes  in  lattice  parameters,  such  as  a  strong 
discontinuous  drop  at  the  lattice  ends,  solitons  were  actually  created  when  the 
parameter  change  was  introduced.  For  example,  in  Figure  V.12,  we  see  a  pair  of 
kinks  just  prior  to  and  just  after  the  introduction  of  a  sinusoidal  modulation  of 
natural  frequency.  When  the  modulation  is  introduced,  a  kink/antikink  pair  arises 
at  element  90  and  leftward  (element  90  is  a  frequency  maximum).  Shortly  after  this, 
the  leftmost  kink  of  the  created  pair  undergoes  an  annihilation  event  with  the 
original  kink,  leaving  behind  a  kink  of  similar  type  located  where  the  original  kink 
would  eventually  have  moved  in  any  case  ~  at  the  frequency  maximum.  Whether 
this  is  in  effect  the  mechanism  for  all  kink  motion  remains  to  be  seen. 

Another  way  one  might  attempt  to  explain  the  difference  between  negative 
energy  kinks,  which  move,  and  positive  energy  kinks  and  domain  walls,  which  don't, 
is  to  assert  that  very  sharp  solitons  are  pinned,  whereas  broad  solitons  are  free  to 
move.  The  ambiguity  arises  from  the  unfortunate  reality  that  only  a  limited  number 
of  cases  have  been  tried  so  far,  and  for  all  of  the  domain  walls  and  positive  energy 
kink  cases,  the  structures  are  very  sharp.  It  may  be  that,  when  a  broad  domain  wall 
is  obtained  (if  they  exist),  it  will  move  exactly  as  the  negative  energy  kinks  do. 
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Fig  V.12.  Transition  by  kink  creation  and  annihilation. 
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Finally,  a  quick  look  at  breathers  revealed  that  they  do  not  seem  to  move  at 
all,  or  even  to  be  affected  by  the  variations,  except  in  this  peak  amplitude.  This  may 
again  be  due  to  sharpness  effects;  certainly  the  envelope  of  a  breather  is  a  rapidly 
changing  function.  This  result  is  striking,  for  breather  motion  has  been  observed 
experimentally  in  water  trough  experiments  (Wu,  et  al.  [1984]). 

The  best  that  can  be  said  after  this  very  preliminary  look  at  the  behavior  of 
nonlinear  lattices  with  nonuniformities  is  that  the  behavior  promises  to  be  as  rich 
and  varied  as  the  uniform  lattice,  with  many  new  phenomena  undoubtedly  awaiting 
later  researchers.  It  would  be  very  interesting  to  study  the  effects  of  nonuniformities 
in  two  dimensional  lattices,  but  that's  another  entire  area  of  research.... 

C.      THE  TODA  LATTICE. 

In  some  ways,  the  most  important  contribution  this  thesis  makes  to  the  study 
of  dynamical  systems  is  the  introduction  of  a  highly  interactive  and  readily  adaptable 
modeling  tool  for  oscillating  systems.  The  modeling  should  be  limited  to  oscillating 
systems,  since  the  errors  incurred  in  using  a  simple  Euler's  method  derivative  may 
accumulate  in  other  types  of  physical  models  (Appendix  A).  The  program  offers  the 
investigator  quick  access  to  an  accurate  model  of  his  system  of  interest  which  allows 
him  to  view  system  dynamics  in  the  time  domain,  frequency  domain,  and  phase 
domain.  Additionally,  a  full  suite  of  file  handling  routines  is  provided  so  the 
investigator  can  save  states  of  his  system,  or  he  can  save  sequences  of  time  sampled 
data  or  spectral  data  in  files  while  observing  system  behavior. 


156 


In  order  to  demonstrate  the  utility  and  flexibility  of  the  program,  a  model  of 
the  well  known  Toda  lattice  (Toda  [1971])  was  developed.  It  required  approximately 
fifteen  minutes  to  completely  alter  the  program,  converting  it  from  a  model  of  one 
physical  system  to  another!  At  the  conclusion  of  that  time,  I  was  taking  data  exactly 
equivalent  to  the  types  of  data  used  throughout  this  thesis.  The  Toda  lattice  is  a 
lattice  where  nearest  neighbors  interact  according  to  an  exponential  interaction 
potential 


*(x)-V*+ax,  V.C.I 

b 


which  yields  an  equation  of  motion  (Kuusela  and  Hietarinta  [1990]) 

x  -o(e~^"~;Vl)-e~*(^r**>).  vc-2 

Here  xn  is  the  usual  variable  representing  departure  from  equilibrium  position  of  the 
nth  element. 

The  free  system  given  by  (V.C.2)  suffers  from  the  same  problem  that  all  free 
systems  do  in  numerical  work  ~  it  is  often  difficult  to  distinguish  solitons  from 
background  radiation,  which  of  course  remains  for  all  time  due  to  the  lack  of 
damping.  Accordingly,  we  desire  to  stabilize  the  structures  of  interest  by  considering 
damped  and  driven  systems.  In  the  case  of  the  Toda  lattice,  there  are  many  possible 
ways  to  add  damping  and  drive.    Geist  and  Lauterborn  [1988]  used  a  sinusoidal 
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driving  force  acting  only  on  the  first  element  of  the  lattice,  and  viscous  damping. 
Since  this  is  quite  distinct  from  the  parametric  drives  used  in  the  lattice  models  we 
have  studied  so  far,  this  method  was  not  used.  Instead,  three  different  drive  schemes 
were  developed.  The  first  two  consisted  simply  of  parametric  drive  via  each  of  the 
two  parameters  of  the  system  (a  and  b);  the  last  was  more  elaborate.  In  this  scheme, 
the  lattice  elements  were  each  placed  in  an  external  potential  well  analogous  to  that 
provided  by  gravity  in  our  previous  lattices;  the  parametric  drive  was  viz  the  natural 
frequency  of  the  element  resulting  solely  from  its  interaction  with  the  potential  well 
(in  other  words,  the  same  parametric  drive  was  used  that  was  used  in  the  other 
models  of  this  thesis). 

Mathematically,  the  first  two  methods  are  given  by 

xn-(a+2r)cos2ut)(e~blx,'~x-l)-e~*x*"~x'))  VC-3 


and 


x  -^-*x"-x--l)-*-*x--x")),  V.C.4 


where 
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B-b+2r\cos2(Dt.  V.C.5 


The  third  method  is  given  by 


-^.-Vi^-Kvi-*^/    *  o«„^,.,rtr  V.C.6 


x-a \e      "   """  -e 


)-(oi0+2r\cos2dit)xn 


As  before,  we  arbitrarily  set  the  natural  frequency  to  1  in  the  program,  without  any 
loss  of  generality. 

All  three  of  these  methods  were  coded  up  and  run  on  the  computer.  The  first 
two  have  proven  to  be  very  difficult  to  control,  with  the  slightest  changes  in  drive 
or  dissipation  producing  large  results.  It  was  possible  to  get  some  stable  states, 
however,  with  one  of  them  shown  in  Figure  V.13.  There  are  two  lattices  shown  in 
the  figure,  because  the  stable  state  is,  more  accurately,  a  quasiperiodic  state  where 
the  two  states  shown  alternate  back  and  forth,  with  a  transient  during  the  alteration 
which  shows  beating  between  the  left  and  right  hand  sides  of  the  lattice.  This  cyclic 
behavior  is  stable,  hence  the  lattice  could  be  termed  quasistable.  However,  it  seems 
to  me  that  the  main  utility  of  the  "pure"  Toda  lattice  is  its  mathematical  tractability 
(it  is  exactly  integrable  and  has  been  studied  extensively  in  the  mathematical 
literature);  it  is  difficult  to  imagine  a  physical  system  which  obeys  the  equations  of 
motion. 

More  physically  interesting  is  the  third  type  of  Toda  lattice  model,  which  can 
be  interpreted  physically  as  particles  with  exponential  interaction  potential  resting 
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Fig.  V.13.  Two  states  of  quasiperiodic  Toda  lattice  with  gamma  drive. 
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each  in  its  own  external  potential  well.  A  physical  realization  of  this  is  the  solid 
state  (this  interpretation  is  due  to  Larraza).  The  exponential  interaction  potential 
is  still  questionable,  but  it  is  at  least  a  feasible  model.  The  third  type  of  model  was 
also  the  most  practical  numerically,  in  the  time  I  was  able  to  allot  to  the  Toda  lattice 
study.  It  should  be  noted  that  it  may  merely  be  a  poor  choice  of  parameters  that 
resulted  in  the  difficulty  in  using  the  first  two  models,  and  a  more  detailed  study  will 
undoubtedly  prove  useful.  The  first  result  obtained  with  the  Toda  lattice  in  an 
external  potential  was  a  negative  energy  kink  similar  to  the  cutoff  kinks  in  the  NLS 
lattice,  in  a  travelling  upper  cutoff  mode.  That  is,  the  mode,  with  its  kink  in  train, 
moved  from  right  to  left  continuously,  at  a  fixed  rate.  This  kink,  could  not  be 
stabilized  as  a  standing  wave  for  system  parameters  close  to  the  original  parameters 
which  yielded  the  results,  although  there  is  no  reason  to  conclude  that  standing  kinks 
might  not  exist  elsewhere  in  the  parameter  space.  In  fact,  for  similar  reasons  we 
cannot  really  conclude  that  there  are  not  travelling  wave  solitons  in  the  nonlinear 
lattice  model  that  the  rest  of  this  thesis  addresses;  we  just  haven't  observed  then  yet. 
Figure  V.14  shows  the  profile  of  this  travelling  wave  domain  wall  at  one  instant  in 
time.   It  maintains  this  form  as  it  moves  to  the  left. 

The  results  shown  in  Figures  V.13  and  V.14  are  only  the  tip  of  the  iceberg,  and 
are  shown  more  to  demonstrate  the  flexibility  of  the  program  and  the  computational 
ease  with  which  a  new  line  of  inquiry  can  be  undertaken.  For  example,  it  was  found 
that  lattices  which  were  driven  in  the  coupling  term  were  very  difficult  to  control 
and  exhibited  at  best  very  limited  stability.  In  any  case,  one  of  the  beauties  of  this 


161 


(0 

^ 

£ 

CO 

c 

c 

0 

•   <M«> 

o 

CO 

Q. 

E 

"CO 

o 

Q 

0 

*- 

X 

U) 

0 

c 
co 

(D 

c 

> 

0 

03 

O 

c 

H 

_co 

co 

CO 

"D 

O 

O 

1- 

H 

i    I 

■ 

■ 
■ 

II 

■ 

■ 

■ 
■ 

■ 
■ 

■ 

■ 

- 

■ 

■ 

■ 
■ 

■ 
■ 


c 
o 

CO 

o 

Q_ 
0 

o 

'■& 

CO 


in      CM 
c\i 


in      t- 


m      o      to      i-      io 

O  O  i- 


CVJ 

I 


If) 
c\i 


(qje)  epnjndujv 


Fig.  V.14.  Travelling  domain  wall  in  Toda  lattice  with  external  potential. 
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work  is  that  the  limiting  factor  in  one's  progress  is  one's  understanding  of  the  physics 
of  lattices;  there  is  no  reason  to  lose  productivity  to  the  computer,  since  the 
modeling  technique  is  well  developed  and  very  easy  to  use  and  mosify.  It  is  to  be 
hoped  that  a  follow  on  thesis  student  will  conduct  a  thorough  study  of  all  three  Toda 
lattice  models  (and  perhaps  others),  in  conjunction  with  an  attempt  to  understand 
the  theory  of  Toda  lattice  solitons.  There  are  some  results  in  this  area  in  the 
literature,  but  very  little  is  known  compared  to  the  simpler  lattice  models. 
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VI.  CONCLUSIONS  AND  PROSPECTS  FOR  FUTURE  WORK 

It  should  be  quite  evident  at  this  point  that  the  behavior  of  nonlinear  lattices 
is  extremely  complex.  We  have  shown  that  kinks  exist  in  virtually  all  modes  of 
lattice  vibration,  and  domain  walls  appear  to  exist  between  all  compatible  modes  (by 
compatibility,  we  mean  all  modes  where  the  higher  amplitude  mode  has  a  negative 
energy  kink  and  the  lower  amplitude  mode  has  a  positive  energy  kink).  Neither  of 
these  two  assertions  has  been  rigorously  proven,  but  the  evidence  for  them  is 
encouraging.  Additionally,  breathers  in  the  cutoff  modes  were  found  to  be  very 
robust,  and  the  remarkable  self -focusing  phenomenon  discovered  experimentally  by 
Denardo  [1990]  was  found  numerically  as  well.  Breathers  in  the  Toda  lattice  and 
in  the  two  dimensional  analog  of  the  basic  lattice  we  discussed  were  also  found  with 
relative  ease,  further  suggesting  that  breathers  are  common  phenomena. 

What  needs  much  more  work,  however,  is  the  theory  underlying  these  many 
nonlinear  structures.  The  NLS  theory  for  the  cutoff  modes  is  satisfactory,  but  it 
covers  only  the  endpoints  of  the  dispersion  curve.  It  is  hoped  that,  underlying  this 
broad  range  of  interrelated  phenomena,  there  may  exist  a  common  equation  of 
evolution  which  will  prove  more  useful  than  the  current  favorite  (the  Korteweg-de 
Vries  Equation),  which  has  limited  physical  application.  If  such  a  general  treatment 
that  includes  the  NLS  as  a  special  case  and  which  treats  all  intermediate  modes  can 
be  found,  it  would  find  broad  application  in  the  physical  sciences,  since  the 
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underlying  equations  of  motion  of  the  systems  considered  here  are  frequently 
encountered.  Examples  include  such  solid  state  topics  as  ferromagnetism, 
superconductivity,  and  nonlinear  optics,  as  well  as  topics  in  plasma  physics, 
cosmology,  and  particle  physics. 

There  are  many  questions  raised  during  the  course  of  this  work  which  can  be 
profitably  addressed  by  follow  on  researchers  using  the  model  developed  here; 
resolution  of  these  questions  will  undoubtedly  ease  the  task  of  achieving  the 
theoretical  understanding  desired.  It  is  still  an  open  question  whether  domain  walls 
and  kinks  are  distinct  fundamental  structures  or  whether  domain  walls  are  always 
matched  kinks  (the  latter  is  my  view,  but  it  is  by  no  means  unanimously  accepted). 
In  addition,  there  are  several  classification  problems  concerning  kinks  in 
intermediate  modes,  which  may  be  useful  to  address  as  a  means  to  guide  the  work 
of  the  theoreticians.  And  it  is  not  yet  known  whether  there  exist  any  breathers  in 
the  intermediate  modes.  Indeed,  the  distinction  between  breathers  and  kinks,  which 
was  made  during  the  study  of  cutoff  mode  solitons  when  the  distinction  was  clear, 
may  not  be  a  useful  distinction  in  the  intermediate  modes.  And,  it  is  unclear  why 
the  NLS  theory  works  so  well  far  outside  the  limits  of  its  validity. 

This  work  focused  mainly  on  kinks  and  domain  walls,  after  the  initial 
verification  of  breathers'  existence  according  to  the  NLS  theory.  It  is  left  to  the  next 
student  (who  has  already  commenced  work)  to  explore  breathers  in  greater  detail; 
it  is  hoped  that  the  results  of  this  later  work  will  dovetail  with  those  presented  here 
to  clear  up  the  basic  questions. 
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We  only  scratched  the  surface  in  this  work  in  the  study  of  two  dimensional 
lattices  and  of  lattices  with  different  equations  of  motion  (i.e.,  the  Toda  lattice).  The 
work  done  was  done  primarily  to  flex  the  muscles  of  the  computer  program  to  verify 
the  ease  of  adaptation  which  was  one  of  the  design  goals  of  the  model.  The  results 
that  were  obtained,  however,  were  very  interesting,  and  make  the  pursuit  of  these 
lines  of  inquiry  highly  desirable.  In  particular,  the  two  dimensional  case  provides  an 
opportunity  to  explore  the  effects  of  diffraction  on  solitons.  Just  from  our 
preliminary  look,  this  effect  appears  to  be  very  important,  as  it  accounts  for  the 
unexpected  stability  of  two  dimensional  breathers  (which  were  thought  before  to  be 
impossible  to  obtain). 

In  conclusion,  we  were  fortunate  in  this  research  effort  to  be  able  to  make 
many  interesting  discoveries  and  to  bring  together  conclusively  some  theoretical  and 
experimental  work.  The  results  obtained  have  great  promise  and  have  opened  up 
several  lines  of  inquiry.  In  addition,  the  value  of  interactive  modeling  was  shown, 
and  a  simple,  fast  model  was  developed  which  can  be  adapted  in  a  matter  of  minutes 
to  model  any  vibrating  discrete  systems.  It  could,  if  desired,  be  used  to  model  wave 
motion  in  continua  as  well,  constituting  in  that  case  a  simple  finite  element  model. 
All  that  is  required  is  that  the  researcher  correctly  enter  the  equations  of  motion  of 
the  system  to  be  studied,  and  to  change  the  I/O  section  if  needed  to  reflect  those 
parameters  only  which  are  relevant  to  the  new  work.  This  moel  may  in  the  end 
prove  to  be  the  most  important  contribution  of  this  research;  in  any  case,  it  is  hoped 
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that  it  will  prove  useful  to  many  future  researchers,  to  whom  it  is  freely  offered  for 
modification  and  use. 
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APPENDIX  A.  NUMERICAL  METHODS. 

Since  all  of  the  phenomena  studied  in  this  thesis  have  been  oscillatory,  it  was 
possible  to  use  a  very  simple  numerical  algorithm,  which  is  a  modification  of  the 
classical  Euler's  method  given  by 

*Jt+Q-*Jto+toJfc  <107) 

and 

\(t+h)'Vn(t)+han(t).  (108) 


where  h  is  the  time  step  used.  Since  the  system  which  I  intended  to  model  was  one 
where  the  rate  of  change  of  x  depended  at  any  given  time  on  the  values  of  x  of  the 
adjacent  elements,  it  seemed  reasonable  to  break  the  calculation  of  (1)  up  as  follows: 

1.  Calculate  the  acceleration  felt  by  each  element,  based  on  the  positions  of 
its  neighbors. 

2.  Update  the  velocities  of  each  element  by  multiplying  the  acceleration  by 
the  time  step. 

3.  After  all  velocities  have  been  updated,  calculate  the  new  positions 
of  the  elements,  and  begin  again  at  1. 
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Mathematically,  this  process  can  be  represented  by  a  pair  of  equations: 

and 

xfi+h)-xJto+k»JLt+h).  (110) 

This  simple  method  turned  out  to  work  very  well  indeed.  It  was  found  that  the 
free,  linear  lattice  behaved  as  expected,  and  that  the  nonlinear  driven  lattice  behaved 
as  hoped.  In  order  to  check  the  quality  of  the  numerics,  an  energy  calculation  was 
included  in  the  program.  The  total  energy  of  a  free  nonlinear  lattice  with  linear 
coupling  between  nearest  neighbors  is  given  by 

It  was  found  in  all  cases  that  the  total  energy  of  the  system  oscillated  about  some 
fixed  value,  with  the  amount  of  oscillation  being  directly  proportional  to  the  time 
increment  used.  Thus  it  was  possible  to  use  a  large  time  step  to  observe  qualitative 
changes  in  system  behavior  more  expeditiously,  and  then  to  switch  to  a  small  step 
size  (typically  0.005T,  where  T  is  the  period  of  the  oscillations)  when  precision  is 
desired  for  quantitative  comparisons  to  theory  or  other  results. 
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As  it  turned  out,  this  method  and  the  conclusions  about  its  efficacy  had  already 
been  developed  (Cromer  [1981]),  so  I  regretfully  have  to  identify  the  method  as  the 
Euler-Cromer  method...  Also,  the  method  solves  a  question  posed  by  Goedde, 
Lichtenberg,  and  Lieberman  [1990],  who  pointed  out  that,  if  too  large  a  time 
discretization  is  used,  parametric  instabilities  develop  in  the  discrete  Sine-Gordon 
Equation  (of  which  our  model  is  a  third  order  approximation).  They  wondered,  in 
their  conclusions,  whether  it  was  possible  to  have  stable  solitons  in  a  discrete  Sine- 
Gordon  system,  or  whether  the  energy  of  the  system  would  always  be 
equipartitioned  among  all  of  the  elements  when  the  time  discretization  was  finite. 
The  results  presented  in  the  body  of  this  thesis,  while  not  a  rigorous  proof,  do 
strongly  suggest  that  it  is  possible  to  have  stable  soliton  solutions  in  a  discrete 
system,  provided  the  time  step  is  reasonably  small  and  a  suitable  method  is  used. 
Then,  the  errors  are  bounded,  and  one  obtains  stable  results  to  at  least  many 
thousands  of  oscillator  periods,  so  that  even  if  the  stability  is  only  asymptotic,  the 
model  is  practically  stable. 

Whereas  unfortunately  the  Euler-Cromer  method  was,  if  somewhat  obscure,  not 
original,  the  vehicle  in  which  it  was  used  is,  I  believe,  original  and  of  great  potential 
utility  to  future  researchers.  In  a  sharp  deviation  from  the  type  of  numerical  analysis 
which  previous  computer  generations  demanded,  wherein  one  started  the  model  with 
some  set  of  initial  conditions  and  parameter  values  and  then  waited  for  the  outcome, 
repeating  the  process  numerous  times  in  order  to  build  up  a  coherent  picture  of 
physical  processes,  this  model  system  is  highly  interactive.  This  allows  the  researcher 
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to  observe  system  behavior  in  "real  time",  and  to  adjust  parameters  as  desired  to 
control  system  behavior  or  to  more  fully  understand  it.  Essentially,  the  model  I  have 
written  is  analogous  to  a  complete  laboratory  setup,  with  data  taking  expedited  and 
unwanted  noise  virtually  eliminated  (except  numerical  "thermal"  noise,  which  is  small 
in  most  cases).  One  can  adjust  the  damping,  drive  amplitude,  or  drive  frequency  via 
either  a  fine  or  coarse  adjust  "knob",  the  knob  being  single  keystrokes  entered  while 
the  model  is  running  at  full  speed.  One  may  also  "kick"  any  chosen  elements  by  any 
desired  amount,  or  "pin"  an  element  at  zero  amplitude,  in  order  to  either  perturb  the 
system  significantly  or  to  move  solitons  along  the  lattice  (as  was  done  in  the 
experimental  studies  with  pendulum  lattices  and  shallow  water  channels). 

The  system's  behavior  can  be  observed  in  the  time  domain  either  by  watching 
the  lattice  on  a  full  screen  display  of  moving  dots,  or  on  a  waterfall  display  where 
trends  are  more  readily  apparent  (but  detail  is  less  clear).  A  zoom  feature  allows 
one  to  rescale  with  the  lattice  dynamics,  or  even  to  measure  numerical  "temperature" 
by  zooming  in  on  nodes  in  modes  such  as  the  "  +  0-0"  mode  (see  Chapter  IE),  until 
the  random  noise  caused  by  roundoff  errors  is  detected.  Several  times  this  was 
done,  and  the  results  were  a  numerical  noise  so  low  that  one  had  to  zoom  240  times 
to  see  it!  A  user  aid  is  available  to  monitor  the  system's  stability,  which  frees  the 
user  to  some  degree  from  having  to  continuously  monitor  the  lattice. 

Additionally,  the  system's  dynamics  can  be  viewed  real  time  in  the  frequency 
and  phase  domains.  In  the  former  case,  a  running  spectrum  (consisting  of 
consecutive  2048  point  Fast  Fourier  Transforms  adapted  from  Press,  Teukolsky  et. 
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al.  [1988])  of  a  single  element  is  displayed  on  the  screen,  along  with  all  system 
parameters.  In  the  latter,  a  phase  portrait  of  up  to  five  elements  at  a  time  is 
displayed,  with  Poincare  section  optional.  In  order  to  eliminate  the  drive  frequency 
and  thus  clean  up  the  display  to  make  interpretation  more  straightforward,  a  rotating 
set  of  axes  was  used,  with  the  coordinate  and  momentum  axes  rotating  synchronous 
with  the  drive  frequency.  In  such  a  display,  a  simple  harmonic  oscillator  appears  as 
a  small  ellipse  (which  vanishes  as  time  increment  approaches  zero). 

The  net  result  of  the  interactive  features  the  model  provides  is  that  one  is  able 
to  develop  intuition  about  system  behavior  as  various  parameters  are  varied,  and  one 
can  experiment  with  various  types  of  perturbations  to  the  system  (including,  as  well 
as  those  just  described,  varying  parameters  randomly,  via  a  gradient,  or  sinusoidally 
(Chapter  IV),  or  just  imposing  a  sudden  ramp  in  amplitude  of  5%  to  test  stability  of 
a  given  configuration).  Such  interaction  leads  one  to  find  structures  which  might 
otherwise  escape  notice,  since  they  often  cannot  be  obtained  simply  by  choosing  a 
set  of  initial  conditions  and  letting  the  system  go.  Additionally,  having  the  variety 
of  representations  available  at  all  times  and  in  real  time  gives  great  diagnostic 
capability  to  the  experienced  researcher. 

Finally,  the  program  was  designed  for  maximum  adaptability  to  the  study  of 
other  problems  concerning  harmonic  systems.  In  Appendix  B,  a  detailed  manual 
describes  not  only  how  to  use  the  program,  but  how  to  alter  it  to  model  other 
systems.    As  discussed  in  Chapter  IV,  a  working  model  of  the  Toda  lattice  was 
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available  in  fifteen  minutes,  with  the  complete  range  of  interactive  tools  available 
for  the  study  of  that  very  challenging  dynamic  system. 
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APPENDIX  B.  PROGRAM  MANUAL. 

A.      INTRODUCTION  AND  GETTING  STARTED. 

The  program  LAl  llCE.C  is  intended  to  be  an  adaptable  modelling  system  for 
oscillating  physical  systems.  The  numerical  method  used,  the  Euler-Cromer  method, 
is  described  in  Appendix  A,  and  should  be  understood  by  the  researcher.  The 
purpose  of  this  appendix  is  to  provide  a  readable  manual  which  describes  how  to  use 
and  modify  the  program  in  order  to  study  physical  systems  using  the  IBM  PC 
(preferably  486  based,  for  speed).  The  program  requires  a  VGA  controller,  and  is 
intended  to  work  in  conjunction  with  screen-dump-to-printer  programs  such  as 
EGALASER.COM,  which  is  commercially  available.  The  first  part  of  this  manual 
describes  the  use  of  the  program,  whereas  the  second  is  addressed  to  those  who  wish 
to  make  alterations  to  it  in  order  to  modify  the  model  being  studied.  Changes  simply 
to  improve  utility  or  to  add  functions  are  not  addressed,  although  there  are  certainly 
many  such  changes  which  might  be  desirable. 

To  start  the  program,  simply  type  LATTICE  at  the  DOS  prompt.  If  you  intend 
to  use  a  data  file  for  initial  conditions,  have  its  name,  including  extension,  ready 
when  you  start  the  program.  The  program  prompts  the  user  for  all  information 
needed  to  start  the  program;  interactive  commands  used  once  the  program  is  started 
are  not  prompted.  However,  as  discussed  in  the  next  section,  the  user  may  be 
prompted  for  additional  information  from  time  to  time.    If  you  desire  to  run  the 
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program  in  conjunction  with  a  screen  dump  program  which  is  TSR  (terminate  and 
stay  resident  -  it  remains  in  your  system's  RAM  after  you  start  it),  you  should  start 
the  TSR  program  before  starting  the  model,  since  you  cannot  do  it  from  within  the 
model. 

B.      USING  THE  PROGRAM. 

After  you  start  the  program,  you  will  be  asked  whether  you  wish  to  use  a  data 
file  for  initial  conditions.  This  is  preferable,  since  otherwise  you  have  little  freedom 
concerning  initial  conditions  (unless  you  modify  the  program  yourself).  If  you  do 
want  to  use  a  data  file,  type  "y"  or  "Y"  and  <  ENTER  >.  The  program  then  asks 
whether  you  want  to  use  the  old  format  data  file  or  not.  The  old  format,  which  have 
been  conventionally  given  the  extension  ".LAT",  consist  of  a  set  of  six  parameters, 
followed  by  two  numbers  for  each  element.  The  first  number  for  each  element  is 
its  amplitude;  the  second  is  its  velocity  (which  is  usually  close  to  zero).  If  this  is  the 
type  of  data  file  you  have,  type  "y".  Otherwise,  your  file  should  be  formatted  with 
only  five  initial  parameters,  but  with  each  element  having  four  numbers.  These 
numbers  represent,  in  order,  the  coupling,  natural  frequency,  amplitude,  and  velocity 
of  the  element.  Typing  any  character  other  than  "y"  at  the  prompt  will  result  in  the 
program  trying  to  read  this  format.  Finally,  the  program  asks  for  the  file  name.  The 
maximum  size  is  30  characters,  so  be  sure  the  file  name,  including  path  if  that  is 
different  from  the  path  your  program  resides  at,  does  not  exceed  thirty  characters! 
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If  you  choose  not  to  use  a  file  for  initial  conditions,  the  program  will  prompt 
you  for  all  of  the  parameters  it  needs.  The  program's  default  initial  conditions  are 
an  amplitude  modulated  cutoff  mode.  If  the  sign  of  the  nonlinear  coefficient  is 
negative,  corresponding  to  a  hardening  lattice,  the  program  sets  up  an  upper  cutoff 
lattice;  otherwise  it  sets  up  a  lower  cutoff  lattice.  The  mode  is  modulated  by  one 
complete  wavelength  of  amplitude  chosen  by  the  user  ("modulation  amplitude",  as 
opposed  to  "mode  amplitude",  which  is  the  amplitude  of  the  cutoff  mode  itself). 

After  either  the  file  data  is  read  or  the  model  parameters  have  been  obtained 
form  the  user,  the  program  asks  you  to  enter  the  time  step.  The  number  you  enter 
will  be  the  time  step,  after  it  is  multiplied  by  one  half  per  cent  of  the  period  of  the 
drive  (this  requires  that  you  do  not  enter  zero  for  drive  frequency,  even  if  you  want 
an  undriven  system,  since  doing  so  will  give  a  divide  by  zero  error,  and  this  program 
doe:  not  include  error  checking).  Typical  choices  for  time  step,  as  entered,  will  be 
from  one  half  to  four,  depending  on  whether  you  are  most  interested  in  exact 
amplitudes  (use  small  time  step)  or  in  getting  relatively  quick,  qualitatively  accurate 
results  (use  large  time  step).  Use  of  time  steps  greater  than  four  does  not  cause  an 
error,  but  it  should  be  done  only  after  the  user  is  familiar  with  the  behavior  of  the 
model  he  is  studying.  If  the  model  is  near  any  sort  of  boundary  in  behavior  type,  or 
if  it  is  near  a  potential  energy  maximum,  using  a  large  time  step  may  perturb  the 
system  into  undesirable  or  even  catastrophic  behavior.  On  the  other  hand,  judicious 
use  of  very  large  time  steps  can  push  the  system  into  stable  behaviors  which  might 


176 


otherwise  have  been  missed  -  experience  is  important  when  choosing  the  time  step! 

1.      The  Default  Graphics  Mode. 

At  this  point,  the  screen  will  go  blank.  Type  ^G  to  get  a  graphics  display 
of  the  system  in  real  time  (~G  means  <  CONTROL  >  G,  and  this  type  of  shorthand 
will  be  used  from  now  on).  The  ^G  command  will  always  do  this,  so  it  is  a  useful 
one  to  memorize  immediately.  If  the  display  is  going  off  the  screen  due  to  excessive 
amplitude,  you  can  "zoom"  out  by  typing  "o".  This  means  in  effect  that  you  step 
back  from  the  system,  so  that  the  observed  maximum  amplitude  increases  (by  a 
factor  of  two,  as  it  happens).  Similarly,  "i"  zooms  you  back  in  by  a  factor  of  two. 
This  is  useful  if  the  amplitudes  get  too  small  to  see  (your  picture  is  limited  by  the 
pixel  resolution  of  the  graphics  mode).  The  graphics  display  will  only  show  40 
elements  on  the  screen;  if  the  lattice  you  are  using  contains  more  than  that,  the  forty 
visible  elements  are  the  middle  elements.  The  other  elements  are  still  being 
calculated  and  will  be  included  if  you  save  your  state  to  a  data  file,  you  simply  can't 
see  them.  The  maximum  number  of  elements  is  150. 

As  a  user  aid,  a  stability  checking  routine  is  included  in  the  program.  This 
routine,  which  is  called  by  typing  "S"  from  any  point  in  the  program,  asks  the  user 
to  select  an  element  for  monitoring  and  then  monitors  that  element  to  see  if  it  is 
"stable".  When  choosing  an  element,  remember  that,  in  the  C  language,  the  first 
element  of  an  array  is  the  zeroth  element.  If  you  want  to  monitor  the  first  element, 
tell  the  computer  to  monitor  element  zero.  Stability  is  determined  by  comparing 
successive  peaks  of  the  amplitude  of  the  chosen  element.   If  twenty  peaks  in  a  row 
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are  all  within  one  percent  of  each  other,  the  system  is  called  stable.  This  does  not 
mean  it  is  stable  in  any  scientific  or  mathematical  sense;  this  function  is  solely  to  aid 
the  user  as  he  sees  fit.  When  this  option  is  chosen,  the  entire  lattice  changes  color 
to  red,  except  the  element  being  monitored  (which  remains  white).  After  stability 
has  been  verified,  the  lattice  returns  to  the  white  color,  except  the  chosen  element, 
which  then  becomes  red. 

The  stability  checker  is  one  of  two  routines  that  uses  the  amplitude  peak 
detection  routine  ~  the  other  is  the  Waterfall  Display,  which  is  described  later.  The 
stability  checker  must  be  active  if  it  is  desired  to  save  the  current  model  state  to  a 
data  file,  since  the  states  are  saved  only  when  the  monitored  element  reaches  an 
amplitude  peak.  If,  while  the  stability  checker  is  active,  you  type  "p",  the  program 
will  pause  and  ask  you  whether  you  wish  to  save  the  state  to  a  file.  If  you  do,  type 
"y"  or  "Y",  and  then  give  the  file  name.  It  is  recommended  that  a  personal  file 
naming  protocol  be  developed  early  so  that  confusion  does  not  occur  when  the 
program  suddenly  asks  for  a  file  name.  If  you  type  "p"  before  the  stability  checker 
is  inactive,  the  pause  routine  will  be  executed  as  soon  as  the  stability  checker  is  next 
activated.  After  dumping  data  to  a  file,  or  declining  to  do  so,  the  program  prompts 
for  a  new  time  step.  The  old  one  is  displayed  for  convenience.  If  you  want  to  keep 
it,  you  must  type  it.  The  program  uses  whatever  number  you  type  as  the  new  time 
step.  Thus,  the  pause  routine  provides  the  only  means  for  changing  the  time  step 
during  model  operation.  The  user  is  cautioned  against  changing  the  time  step  while 
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the  Fourier  Analysis  routine  is  operating  (see  below),  since  a  change  in  the  time  step 
then  will  corrupt  the  Fourier  data  and  give  erroneous  spectra. 

The  drive  parameters  and  the  dissipation  parameter  can  be  changed  while 
the  program  is  operating  (it  does  not  need  to  be  in  Default  Graphics  Mode  to  do 
this).  The  letters  U,V,D,  and  E  represent  coarse  and  fine  increases  and  decreases 
in  the  applicable  parameter.  A  <  SHIFT  >  operation  with  one  of  these  letters  means 
the  parameter  affected  is  drive  frequency.  Thus  U  means  coarse  increase  in  drive 
frequency.  Similarly,  ^U  means  coarse  increase  in  drive  amplitude,  and  u  means 
coarse  increase  in  dissipation.  The  commands  ^ZjZ,  and  z  mean  that  the  indicated 
parameter  (drive  amplitude,  frequency,  and  dissipation,  respectively)  gets  set  directly 
to  zero.  This  is  very  useful  if  your  model  starts  to  blow  up  ~  type  ~Z  and  eliminate 
drive,  so  that  system  energy  will  at  best  remain  constant  (it  will  decrease  if  there  is 
any  damping  still  present).  Finally,  s  brings  all  elements  immediately  to  rest  at  zero 
amplitude  (another  useful,  but  destructive,  way  to  stop  a  system  from  blowing  up). 

It  is  possible  to  interact  on  an  element  by  element  basis  with  your  model 
system,  using  the  ~K  and  ~L  keys.  The  first  of  these  "kicks"  an  element,  which 
means  a  step  change  in  amplitude  is  made.  Velocity  does  not  change,  however, 
which  should  be  kept  in  mind  by  the  observer.  The  second  "pins"  an  element  at  zero 
amplitude  and  velocity.  The  element  stays  pinned  until  *L  is  typed  again  and  the 
same  element  selected  (i.e.,  ~L  is  a  toggle  switch,  as  are  several  other  commands). 

The  lattice  can  be  changed  on  a  global  scale  by  introducing 
nonuniformities  into  the  parameters  representing  coupling  and  natural  frequency, 
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using  the  W  keystroke.  These  parameters  are  stored  as  arrays,  with  a  separate  value 
for  each  element.  This  is  so  they  may  be  made  to  vary  from  element  to  element. 
There  are  three  types  of  variation  available,  and  they  can  be  applied  successively  as 
many  times  as  desired  (this  includes  "undoing"  a  variation  by  introducing  an  identical 
but  opposite  polarity  variation  ~  except,  of  course,  for  random  variations).  The  first 
type  is  a  random,  five  percent  (peak  to  peak)  variation  of  the  parameter  (only  one 
of  the  two  parameters  can  be  varied  for  each  keystroke,  although  both  can  be  varied 
by  successive  keystrokes).  The  second  type  is  a  linear  gradient.  This  gradient  is 
positive,  with  the  zeroth  element's  parameter  remaining  unchanged,  and  the  last 
element's  parameter  being  changed  by  the  full  amount  of  the  gradient.  This  amount 
(in  percent)  is  chosen  by  the  user.  The  third  type  is  a  sinusoidal  variation,  in  which 
the  amplitude  of  the  sinusoidal  variation  and  the  (integer)  number  of  wavelengths 
of  the  modulation  along  the  entire  length  of  the  lattice  are  input  by  the  user.  These 
parameter  variations  are  useful  for  studying  the  effects  of  nonuniformities  on  the 
dynamics  of  the  lattice  being  studied.  Additional  types  of  variations  can  easily  be 
added  by  the  user  (see  modification  instructions  below). 

Another  useful  function,  dumping  a  single  element's  behavior  to  a  data  file 
in  the  form  of  a  pair  of  time  series  representing  displacement  and  velocity,  is 
available  to  the  user  in  all  modes  ("modes"  refers  to  Default  Graphics,  Phase  Plot, 
Text,  Waterfall,  and  Spectrum  display  modes).  Type  ~F  to  start  file  dumping.  The 
program  will  prompt  for  the  element  to  be  dumped;  then  the  program  returns  to 
normal  operation.  To  stop  the  dump,  type  ^F  again,  and  give  the  name  of  the  data 
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file  you  want  to  output  to.  If  you  don't  type  ^F  again,  the  program  will  prompt  you 
after  the  4096  point  buffer  is  full. 

Finally,  typing  n  perturbs  the  system  by  kicking  each  element  by  up  to 
10%,  depending  on  its  location.  This  command  is  useful  when  the  stability  to 
perturbation  of  a  state  is  to  be  determined.  If  the  state  returns  to  normal  after  a 
few  of  these  perturbations,  it  is  definitely  stable! 

2.      Phase  Plot  Mode. 

Phase  plane  diagrams,  or  phase  plots,  are  very  useful  tools  in  gaining  an 
understanding  of  dynamical  systems.  Accordingly,  this  type  of  display  is  available 
to  the  user  via  the  P  command.  This  command,  which  can  be  used  anytime  the 
model  is  running,  shifts  the  display  to  phase  plot  mode.  The  program  asks  whether 
Poincare  sections  are  desired.  Typically,  they  are,  since  they  provide  a  much  cleaner 
and  easier  way  to  interpret  display.  Type  y  to  get  Poincare  sections,  anything  else 
to  get  phase  plots  with  continuous  updates  (one  dot  per  time  step).  The  program 
then  asks  for  up  to  five  elements  to  display.  Entering  five  will  result  in  the  model 
resuming  and  the  phase  plot  beginning  to  build  up.  Entering  fewer  than  five  results 
in  nothing  happening  until  the  user  enters  999  to  indicate  that  no  further  elements 
are  desired  (this  entry  is  prompted,  so  you  don't  have  to  remember  it).  When  in 
phase  plot  mode,  all  of  the  commands  which  affect  the  lattice  are  valid,  and  their 
effects  will  be  seen  in  the  phase  plots.  It  is  recommended  that  you  experiment  with 
this  a  lot,  because  it  is  very  interesting  and  useful,  and  phase  plot  information  is 
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file  you  want  to  output  to.  If  you  don't  type  T  again,  the  program  will  prompt  you 
after  the  4096  point  buffer  is  full. 

Finally,  typing  n  perturbs  the  system  by  kicking  each  element  by  up  to 
10%,  depending  on  its  location.  This  command  is  useful  when  the  stability  to 
perturbation  of  a  state  is  to  be  determined.  If  the  state  returns  to  normal  after  a 
few  of  these  perturbations,  it  is  definitely  stable! 

2.      Phase  Plot  Mode. 

Phase  plane  diagrams,  or  phase  plots,  are  very  useful  tools  in  gaining  an 
understanding  of  dynamical  systems.  Accordingly,  this  type  of  display  is  available 
to  the  user  via  the  P  command.  This  command,  which  can  be  used  anytime  the 
model  is  running,  shifts  the  display  to  phase  plot  mode.  The  program  asks  whether 
Poincare  sections  are  desired.  Typically,  they  are,  since  they  provide  a  much  cleaner 
and  easier  way  to  interpret  display.  Type  y  to  get  Poincare  sections,  anything  else 
to  get  phase  plots  with  continuous  updates  (one  dot  per  time  step).  The  program 
then  asks  for  up  to  five  elements  to  display.  Entering  five  will  result  in  the  model 
resuming  and  the  phase  plot  beginning  to  build  up.  Entering  fewer  than  five  results 
in  nothing  happening  until  the  user  enters  999  to  indicate  that  no  further  elements 
are  desired  (this  entry  is  prompted,  so  you  don't  have  to  remember  it).  When  in 
phase  plot  mode,  all  of  the  commands  which  affect  the  lattice  are  valid,  and  their 
effects  will  be  seen  in  the  phase  plots.  It  is  recommended  that  you  experiment  with 
this  a  lot,  because  it  is  very  interesting  and  useful,  and  phase  plot  information  is 
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most  useful  only  after  significant  experience  has  been  built  up  looking  at  them  under 
various  conditions. 

3.  Text  Mode. 

The  text  mode  is  activated  by  T.  If  this  is  preceded  at  any  time  by  ~H, 
the  text  mode  display  will  conclude  with  a  display  of  total  system  energy,  calculated 
as  described  in  Appendix  A.  In  this  mode,  the  text  mode  functions  as  a  mode  in 
which  model  operation  continues;  in  any  other  case,  the  model  operation  is 
suspended  until  you  leave  text  mode  by  using  the  appropriate  command  to  enter 
another  mode.  The  text  mode  gives  the  values  of  all  system  parameters,  and  then 
waits  for  you  to  hit  a  key.  When  you  do,  it  displays  the  amplitude  and  velocity  of 
each  element.  In  all  honesty,  this  is  not  a  very  refined  mode,  so  it  is  really  useful 
only  for  looking  at  system  parameters  and  system  energy. 

4.  Spectrum  Mode. 

By  typing  ~X,  the  user  enters  Spectrum  Mode.  The  user  is  prompted  to 
enter  the  number  of  the  element  whose  spectrum  is  desired.  Then  execution  returns 
to  the  mode  the  user  was  in  when  he  typed  ~X.  However,  the  amplitude  of  the 
chosen  element  is  now  dumped  every  tenth  time  increment  to  a  data  array  of  2048 
elements.  When  the  array  is  full,  the  Fast  Fourier  Transform  is  calculated  and 
displayed.  This  display  then  remains  on  the  screen  until  it  is  updated  (the  model 
continues  to  run,  and  sequential  FFTs  are  calculated  every  20480  time  increments 
until  the  user  enters  a  different  mode  using  the  appropriate  command).  The  system 
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parameters  and  maximum  frequency  shown  are  displayed  at  the  top  of  the  screen, 
and  a  frequency  axis  is  displayed  at  the  bottom  so  that  numerical  values  of  the 
spectral  components  can  be  estimated.  If  it  is  desired  for  the  complex  FFT  values 
to  be  dumped  to  a  data  file  for  later  analysis,  type  y  and  then  enter  the  filename 
when  prompted. 

5.     Waterfall  Mode. 

The  waterfall  mode,  entered  via  the  ^W  command,  is  useful  for 
monitoring  trends,  particularly  after  some  sort  of  change  has  been  made.  I  found 
it  most  useful  when  I  introduced  lattice  nonuniformities.  If  I  immediately  typed  AW, 
I  got  a  good  visual  indication  of  the  exact  response  of  the  lattice  to  the  step  changes 
I  introduced.  This  type  of  use  is  commended  to  the  user.  When  in  Waterfall  Mode, 
a  waterfall  display  of  110  iterations  is  displayed,  one  iteration  at  a  time.  It  is  best 
to  enter  Waterfall  Mode  with  the  stability  checker  inactive,  to  avoid  multiple 
waterfall  entries  for  the  same  lattice  state.  If  this  is  done,  the  lattice  is  displayed 
once  during  each  cycle,  at  the  same  point  in  the  cycle,  so  that  useful  comparisons  can 
be  made.  When  the  Waterfall  Display  is  full,  the  system  does  not  continue  to 
operate.  This  is  so  that  one  can  pause  to  examine  the  display,  dump  it  to  the  printer 
if  desired,  and  then  press  ~W  again  to  restart  the  model  and  refresh  the  display. 
Operating  in  this  way  makes  it  possible  to  have  a  continuous  series  of  waterfall 
displays  recording  system  behavior  for  long  periods  of  time. 

The  zoom  in  and  out  commands  (i  and  o)  are  operative  in  this  mode  (and 
the  phase  plot  mode),  so  the  display  can  be  adjusted  for  maximum  effectiveness. 
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Also,  typing  ~W  while  the  waterfall  display  is  not  full  toggles  the  display  off.  The 
model  then  continues  to  run,  but  no  new  display  will  appear  until  the  appropriate 
mode  command  (~G,  *T,  P,  ~X,  or  ~W)  is  entered. 

When  finished  with  the  program,  type  ^Q  to  quit. 

C.      MODIFYING  THE  PROGRAM. 

The  program  was  designed  to  be  as  modular  as  possible,  although  there  was  no 
great  effort  expended  to  optimize  performance  in  general,  as  the  emphasis  was  on 
research  and  not  programming  finesse.  This  modularity  makes  it  easy  to  modify  the 
program  to  model  other  physical  systems.  Other  modifications,  such  as  to  add 
additional  features  or  display  modes,  will  not  be  addressed  here. 

To  model  a  physical  system,  the  first  requirement  is  to  develop  the  discrete 
equations  of  motion.  Discreteness  is  required  for  computational  purposes,  although 
time  can  be  treated  as  continuous  (it  is  discrete,  but  the  equations  do  not  need  to 
reflect  that,  since  we  calculate  the  acceleration  and  velocity  at  discrete  instants  in 
time,  thus  implicitly  discretizing  time  for  you).  These  equations  should  then  be 
coded  up  in  the  C  language,  and  then  inserted  directly  in  a  copy  of  the  waterfall 
program  in  place  of  the  equations  of  motion  currently  installed.  The  equation 
entered  should  be  of  the  form 
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xn  -/(coordinates, momenta, t) .  (112) 

As  long  as  this  form  is  used,  and  as  long  as  the  system  described  by  the  equations 
of  motion  is  oscillatory,  at  least  in  steady  state,  then  the  numerical  methods  used 
should  provide  excellent  results.  The  user  should  verify  this  in  every  case  by 
checking  the  numerics  as  described  in  Appendix  A. 

Once  the  equations  of  motion  are  modified,  you  need  to  initiate  any  variables 
or  parameters  that  appear  in  the  equations  of  motion  that  are  not  already  initiated. 
For  example,  when  I  modified  the  program  to  create  a  Toda  Lattice  model  (see 
Chapter  IV),  I  had  to  add  a  variable  b,  representing  the  exponential  coupling 
parameter.  If  these  parameters  are  to  be  input  by  the  user,  then  the  user_init 
function  should  be  modified;  also,  any  references  to  parameters  which  are  removed 
in  the  new  model  should  be  deleted  from  the  I/O  statements  to  minimize  confusion. 
The  parameter  changes  should  also  be  reflected  in  the  file  I/O  sequences  in 
user_init  and  process_pause  functions.  Also,  for  completeness,  the  user  should 
modify  the  parameters  displayed  in  displaytext,  init_phasplot,  and  plotspectrum 
functions.  When  these  changes  are  made,  the  model  is  ready  to  compile  using  either 
the  QUICKC  compiler  or,  if  time  optimization  is  required,  the  Microsoft  C  version 
6.0  compiler.  That  is  all  there  is  to  modelling  a  new  physical  system!  As  I  stated 
in  Chapter  IV,  this  is  a  real  strength  of  this  program;  I  created  a  Toda  lattice  model 
in  fifteen  minutes. 
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SUMMARY  OF  COMMANDS 

O  Quit  the  Program. 

G  Enter  Default  Graphics  Mode. 

T  Enter  Text  Mode. 

~W  Enter  Waterfall  Mode. 

P  Enter  Phase  Plot  Mode. 

~X  Enter  Spectrum  Mode. 

p  Pause  program  (file  I/O,  change  time  step). 

W  Introduce  nonuniformities  into  lattice. 

~H  Monitor  total  system  energy  (in  text  mode  only). 

~F  Toggle  time  series  data  dump  to  file. 

i/o  Zoom  in/out  by  a  factor  of  two  in  amplitude. 

^L  Pin  an  element. 

~K  Kick  an  element. 

y  Dump  spectrum  to  data  file. 

n  Perturn  the  system  by  a  ten  percent  amplitude  gradient. 

S  Toggle  stability  checker  on/off. 

s  Stop  all  elements. 

^O/^I  Increase/decrease  Poincare  strobe  frequency. 

^U^V^D^E  Coarse/fine  increase/decrease  of  drive  amplitude. 

U,V,D,E  Coarse/fine  increase/decrease  of  drive  frequency. 

u,v,d,e  Coarse/fine  increase/decrease  of  dissipation. 
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~Z,Z,z  Zero  out  the  drive  amplitude/frequency/dissipation. 

~R  Reset  and  restart  the  program. 
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C.  PROGRAM  CODE 

/*        PROGRAM  WATERFALL  GENERALIZED  NONLINEAR  LATTICE  MODEL 

VERSION  3.0 

WRITTEN  BY  BRIAN  GALVIN 

LAST_UPDATE  10  JAN  1991 

THIS  PROGRAM  SIMULATES  A  GENERAL  LATTICE  WITH  EQUATIONS  OF  MOTION 
THAT  CAN  BE  SUBSTITUTED  IN  WHERE  INDICATED.   VARIOUS  INTERACTIVE 
FEATURES  ARE  PROVIDED,  WHICH  ARE  EXPLAINED  IN  THE  PROGRAMMER'S 
MANUAL.  A  WATERFALL  DISPLAY  IS  ADDED  TO  MONITOR  SOLITON  MOTION 
DUE  TO  MEDIUM  NONUNIFORMITY  EFFECTS.   COUPLING  AND  NATURAL 
FREQUENCY  CAN  BE  MADE  NONUNIFORM  VIA  SHIFT-W.   */ 

/*  ***  PREPROCESSOR  DIRECTIVES  ***  */ 

linclude  "\quickc\inc\BIOS.H" 

/include  "\quickc\inc\graph.h" 

linclude  "\quickc\inc\stdlib.h" 

linclude  "\quickc\inc\CONIO.H" 

linclude  "\quickc\inc\MATH.H" 

linclude  "\quickc\inc\STDIO.H" 

linclude  "XquickcXincXtime.h" 

/*  ***  MACRO  DEFINITIONS    ***  */ 

Idefine  SQR(a)  ((a)*(a)) 

Idefine  CUB(a)  ( (a  )  * ( a  )  * ( a  )  ) 

Idefine  SWAP(a,b)  tempr=(a) ; (a)=(b) ; (b)=tenpr  /*  USED  IN  FFT  ROUTINE  */ 

Idefine  DOFOR(i.to)  f or( i=0; i<to; i++)  /*  SIMPLIFIED  DO  LOOP  COMMAND  */ 

Idefine  PI  3.14159265359 

Idefine  SCREEN_CORRECTION_FACTOR  1.05 

Idefine  eta_increment  0.01  /*  These  increments  are  used  in  */ 

Idefine  beta_increment  0.01  /*    adjusting  parameters  */ 

Idefine  oraega_increroent  0.001 

Idefine  STABILITY_INCREMENT  0.01 

Idefine  DISPLAY_INCREMENT  DINC 

Idefine  PHASE_INCREMENT  PINC 

Idefine  MIDDLE_ELEMENT  (no_pendulums/2 ) 

Idefine  text_flag  flags[0]  /*  flags  have  names  in  body  of  program,  */ 

Idefine  graphics_f lag  flags[l]      /*  but  they  all  form  a  single  array.   */ 

Idefine  f ile_dump_f lag  flags[2] 

Idefine  stabil ity_f lag  flags[3] 

Idefine  peak_flag  flags[4] 

Idefine  stable_flag  flags[5] 

Idefine  phase_flag  flags[6] 

Idefine  pause_flag  flags[7] 

Idefine  pause_flag2  flags[8] 

Idefine  stop_flag  flags[9] 

Idefine  spectrum_f lag  flags(10] 

Idefine  dump_spectrum_f lag  flags[ll] 

Idefine  energy_flag  flags[l2] 

Idefine  waterf all_f lag  flags[13] 

Idefine  wait.flag  flags[14] 

/*   ***   DYNAMICAL  VARIABLE  DECLARATION.   THESE  ARE  GLOBAL  VARIABLES  ***  */ 

double  coordinate [ 150] , momentum [ 150] ,old_coordinate[ 150) ,old_momentum[  150] ; 

double  accel era t ion, gamma [ 150 ] ,mean_garoraa ,beta,omega0[ 150 ] ,mean_omegaO , 
eta , omega , alpha , model_time, time_int , 

max_amp,mode_amp,time_series_disp( 8000 ] ,phase_elements[5 ] , 
peak_record [ 20 ) , pinned_elements [ 150 ), DINC , PINC , period , 


189 


spectrum [4098] ,temp2[4096] .energy , gradient ,wat_inc; 
int    no_pendulums,f lags[ 15] , counter 1 , phase_counter ,phase_int , f irst_element, 
chosen_element,stability_eleraent,disp_co lor , colors [ 5 ] , counter 2 , 
spectrum_element , spectrum_counter , sample_counter , energy_counter , 
waterf all_counter , color_counter ; 

nain()  ( 

FILE  *fp; 

int  c,i , j,k,l ,m,n,xx; 

double  modulation; 

char  ans [ 5 ] ; 

_clearscreen(_GCLEARSCREEN) ; 


user_init( ) ; 

_setvideomode(_MRES4C0L0R) ; 

_setcolor(l) ; 

xx=0; 

wat_inc=l; 

disp_color=l ; 


/*  INITIALIZE  COUNTERS  HERE  */ 

counter2=0; 
energy_counter=0 ; 
color_counter=0 ; 

/*  ***  BODY  OF  PROGRAM  STARTS  HERE  ***  */ 


/*  STARTS  PROGRAM  WITH  USER  INPUT  */ 


while (stop_flag==0)  ( 
energy=0 ; 
if (pause_flag2==l)  { 

pause_f lag2=0 ; 

pause_f lag=0; 

process_pause( ) ; 
)   /*  if (pause. . . )  */ 
if (kbhit()!=0)  ( 

c=getch( ) ; 

if(c=17)  ( 

stop_f lng-1 ; 

) 

else  if(c==112)  ( 

pause_f lag=l; 

) 

else  if(c==20)   ( 


/*  A  peak  has  been  detected,  so  the  */ 
/*  pause  routine  is  called  */ 


/* 

/* 


CHECK  FOR  USER  INPUT  DURING  MODEL  RUN  */ 
*Q 


/•  P 


/*  *T 


_clearscreen(_GCLEARSCREEH) ; 

text_flag=l; 

graphics_f lag=0; 

waterf all_flag=0 ; 

wait_f lag=0; 

spectrum_f lag=0; 

display_text( ) ; 
) 
else  if(c==7)   (         /*  *G 

_clearscreen(_GCLEARSCREEH) ; 

_setvideomode(_MRES4COLOR) ; 

_setcolor( 1) ; 

phase_f lag=0; 

wait_f lag=0; 

spectrum_f lag=0 ; 

waterf al l_f lng-0; 

text_f lag=0; 


Quit  program  */ 

Pause  program/dump  state  */ 

Text  Mode  */ 


Graphics  Mode  */ 
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graphics_f lag=l ; 
) 
else  if(c==8)   (         /*  AH       :   Monitor  energy  in  text  node  */ 

if (energy_f lag==0)  energy_f lag=l ; 

else  energy_f lag=0; 
) 

/*  •**  THE  FOLLOWING  CASE  IS  FOR  CHANGING  THE  UNIFORMITY  OF  THE 
LATTICE.  COUPLING  OR  NATURAL  FREQUENCY  CAN  BE  VARIED. 
THE  THREE  TYPES  OF  VARIATION  ARE: 

1.  RANDOM  VARIATION  OF  FIVE  PER  CENT 

2.  LINEAR  GRADIENT  OF  ANY  AMOUNT 

3.  SINUSOIDAL  GRADIENT  OF  ANY  AMOUNT,  WITH  AN 

ARBITRARY  NUMBER  OF  FULL  WAVELENGTHS  OVER 
THE  LENGTH  OF  THE  LATTICE   ***  */ 

else  if(c=-87)   {       /*  w     :   Kick  in  parameter  variations  */ 
_setvideomode(_DEFAULTMODE) ; 
srand( (unsigned)time(NULL) ) ; 

printf ("\nEnter  1  if  you  wish  to  vary  coupling:   " ) ; 
scanf ("%d",&j); 
if(j=«l)  { 

printf ("\nEnter  1  (random),  2  (gradient),  or  3(sine)  desired:   "); 
scanf ("%d" ,4k); 
if(k==2)  ( 
printf ( "\nEnter  gradient  (percentage  over  entire  lattice):   "); 
scanf ("%lf ",&gradient) ; 
DOFOR(i ,no_pendulums)  ( 

gamma[ i )=gamma( i ]+mean_garama*i*gradient/( 100*no_pendulums) ; 
)  /*  DOFOR  */ 
)    /*   if(j==2)  */ 
else  if(k==3)  ( 

printf ( "\nEnter  modulation  amplitude  (percentage):  "); 

scanf  ("tlf.&gradient)  ; 

printf ( "\nFnter  integer  number  of  wavelengths  to  use:   "); 

scanf ("\d",Sl ); 

DOFOR ( i ,no_pendulums)   ( 

gamma ( i  J=gamma[ i ) +mean_gamma* gradient/1 00* (-cos ( 2*l*PI*i/ 
no_pendulums ) ) ; 
) 
) 

else  DOFOR( i ,no_pendulums)  ( 
j=rand(  )  ; 

gamma ( i ]=mean_gamma* ( .95+ . 1*  j/RAND_MAX) ; 
) 
) 
else  ( 

printf ("\nVarying  omegao. . .\n" ) ; 

printf ( "\nEnter  1  (random),  2  (gradient),  or  3(sine)  desired:   ") 
scanf ("%d",&k) ; 
if(k==2;  ( 
printf ( "\nEnter  gradient  (percentage  over  entire  lattice):   "); 
scanf ("%lf",6gradient); 
DOFOR( i ,no_pendulums)  ( 

omega0[ i  J=omega0( i ]+mean_omegaO*i*gradient/( 100*no_pendulums) ; 
) 
) 

else  if(k==3)  ( 
printf ( "\nEnter  modulation  amplitude  (percentage):   "); 
scanf ("%lf",&gradient ); 
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printf ( "\nEnter  integer  number  of  wavelengths  to  use:   " ) ; 
scanf ("%dM,&l) ; 
DOFOR( i ,no_pendulums)  { 
omegaO[  i  ]=omega0[i]+inean_omega0*gradient/100*  ( -cos(  2*l*PI*i/no_pendulunis) ) ; 
) 
) 
else  { 

DOFOR(i , no_pendulums)  ( 
j=rand( ) ; 

printf ("Random  number  is:  %d   RAND_MAX  is  %d\n" , j ,RAND_MAX) ; 
omegaO [ i ]=mean_omegaO* ( . 95+ . 1*  j/RAND_MAX ) ; 
) 
> 
) 
) 

else  if(c==23)   (        /*  *W      :   Waterfall  Display  Mode  */ 
if (waterfall_flag==0)  ( 
phase_f lag=0; 
graphics_f lag=o ; 
text_f lag=0; 
spectrum_f lag=0 ; 
waterf all_f lag=l ; 
init_waterf all ( ) ; 

) 
else  ( 

if (wait_f lag==l)  ( 
wait_f lag=0; 
init_waterfall( ) ; 
) 

else  waterf all_flag=0; 
} 
) 

else  if(c==80)   (       /*  P        :   Phase  Plot  Mode  */ 
if (phase_f lag==0)  init_phasplot( ) ; 
else  ( 

phase_f lag=0; 
spectrum_f lag=0; 
graphics_f lag=l  ; 
waterf all_flag=0 ; 
wait_f lag=0; 

_setvideomode(_MRES4C0L0R) ; 
_clearscreen(_GCLEARSCREEH) ; 
)  /*  else  */ 

) 

else  if(c==6)   (         /*  "F        :   Time  series  to  file  */ 
if ( f ile_dump_f lag==0)  start_f ile_dump( ) ; 
else  ( 

stop_f ile_dump( ) ; 
f ile_dump_f lag=0; 
) 
) 

else  if(c==24)   (       /*  *X       :   Calculate  spectrum  */ 
spectrum_f lag=l ; 
_setvideomode(_DEFAULTMODE) ; 
_clearscreen(_GCLEARSCREEN) ; 

printf ("\nEnter  number  of  element  to  be  analyzed:   "); 
scanf ( '^d" , &spectrum_element ) ; 
_setvideomode(_MRES4C0LOR) ; 
_clearscreen(_GCLEARSCREEN) ; 
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else  if(c==21)   (      '  /*  AU 

eta=eta+eta_increment ; 

stability_restart( ) ; 
) 
else  if(c«4)   (        /*  AD 

eta«»eta-eta_increment ; 

stability_restart( ) ; 


else 


else 


) 


if(c=22)   (        /*  AV 
eta=eta+eta_increment/10 ; 
stability_restart( ) ; 

if(c==5)   (        /*  'E 
eta=eta-eta_increment/10 ; 
stability_restart( ) ; 


( 


/< 


else  if(c==26) 

eta=0 ; 

stability_restart( ) ; 
) 
else  if(c==85)   (       /*  u        : 

omega=omega+omega_increment ; 

stability_restart( ) ; 
) 
else  if(c==68)   {       /*  D 

time_int=time_int*200/period; 

omega=omega-omega_increinent ; 

period=2*PI/omega; 

t ime_int=time_int* period/200 ; 

stability_restart( ) ; 
» 
else  if(c==86)   (        /*  V        : 

time_int=time_int*200/period; 

omega=omega+omega_increment/l 0 : 

period=2*PI/omega ; 

time_int=tinie_int*period/200; 

stability_restart( ) ; 
) 
else  if(c==69)   (        /*  E 

time_int=time_int*200/period; 

omega=omcga-omega_increment ; 

per iod=2* PI /omega; 

time_int=time_int*period/200; 

stability_restart( ) ; 


) 
else 


) 
else 

) 


Increase  drive  (coarse)  */ 


Decrease  drive   (coarse)  */ 


Increase  drive   (fine)   */ 


Decrease  drive    (fine)  */ 


Turn  off  drive  */ 


Increase  freguency  (coarse)*/ 


Decrease  freguency  (coarse)*/ 


Increase  freguency  (fine)*/ 


Decrease  freguency  (fine)*/ 


if(c==90)   (       /*  S        : 
omega=0 ; 
stability_restart( ) ; 

if(c==117)   {       /*  u        : 
beta=beta+beta_increment ; 
stability_restart( ) ; 

else  if(c==100) 
beta=beta-beta_increment; 
stability_restart( ) ; 


Set  freguency  to  zero  */ 


Increase  damping  (coarse)*/ 


) 

else  if(c==118)   (       /*  v 

beta=beta+beta_increment/10 ; 

stability_restart( ) ; 


/*  d 


Increase  damping  (fine)*/ 


Deer 
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else  if(c==121)  (        /*  y  :   Dump  spectrum   */ 

if (spectrum_f lag==l )  dump_spectrum_f lag=l ; 


) 

else  if(c-=101)   {      /*  e 

beta=beta-beta_increraent ; 

stability_restart( ) ; 


} 


( 


else  if(c«=l22) 
beta-O; 
stability_restart( ) ; 


/*  2 


) 


/* 


else  if(c==12)   ( 

pin_elements( ) ; 
stability_restart( ) ; 
) 

else  if(c==ll)   {       /*  AK 
kick_elements( ) ; 
stability_restart( ) ; 
) 

else  if(c==105)    ( 
DINC=DINC/2; 
PINOPINC*2; 
wat_inc=wat_inc*2 ; 


/*  i 


Decrease  damping  (fine)*/ 
Turn  off  damping  */ 
Pin  elements  */ 
Kick  elements  */ 
Zoom  in   */ 


) 


/*  o 


Zoom  out 


'/ 


/*   n 


else  if(c==lll)  ( 
DINC=DINC*2; 
PINC=PINC/2; 
wat_inc=wat_inc/2 ; 
) 
else  if(c==HO)  ( 

DOFOR( i ,no_pendulums)  ( 

modulation=. l*i/no_pendulums; 

coordinate! i ]=coordinate[ i ]*(1 .O0+modulation) ; 
) 
) 
else  if(c==15)   (       /*  *o 

phase_int=phase_int-l ; 
) 
else  if(c==9)    {       /*  'I 

phase_int=phase_int+l ; 
) 
else  if(c==l8)   <       /*  *R 

_setvideomode(_DEFAULTMODE) ; 
_clearscreen(_GCLEARSCREEN) ; 
user_init( ) ; 
stability_restart( ) ; 

if (phase_f lag==0)  _setvideomode(_MRES4C0L0R) ; 
else  _setvideomode(_VRES16C0L0R) ; 
) 

else  if(c==ll5)  {  /*  s        :   Stop  everything  as  is  */ 

DOFOR(i ,no_pendulums)  ( 
coordinate! i ]=0; 
momentum [ i ] =0 ; 
)  /*  DO FOR  */ 
stability_restart( ) ; 

> 

else  if(c==83)  (         /*  S         :   Monitor  stability  */ 

monitor_stability( ) ; 
) 

/*  if(kbhit...)  */ 


Perturb  system  */ 


Decrease  strobe  freq  */ 
Increase  strobe  freq  */ 
Restart  program  */ 
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spectrum_counter++ ; 

sample_counter=0 ; 
} 
else  ( 

compute_spectrum() ;   /*  COMPUTE  SPECTRUM  IF  SPECTRUM  DATA  FULL  */ 

_clearscreen(_GCLEARSCREEN) ; 

_setvideoraode(_VRES16COLOR) ; 

plot_spectrura( ) ; 

phase_f lag=0; 

text_f lag=0; 

graphics_f lag=0; 

DOFOR( j,4096)  spectrum[ j J=0 ; 

spectrum_counter=0 ; 

sample_counter=0 ; 
)  /*  else  */ 
}  /*  if  ((spec...  */ 

/*  CHECK  FOR  PEAKS  AND  MONITOR  STABILITY/SET  PAUSE  FLAG/INITIATE 
WATERFAL  DISPLAY  AS  REQUESTED  BY  THE  USER  */ 

if  ( ( ( i=stability_element )  &&  ( stability_f  lag==l )  &&  ( stable_f  lag==0 ) )  | 
( (i=(no_pendulums-l)  )&&(waterf  all_f  lag==l ) ) )  ( 
if ( ( coordinate [ i ]>0)&&( coordinate [i ]<old_coordinate[ i  J ) 
4&(peak_flag==0) )  ( 
peak_f lag=l; 

if (pause_f lag==l )  pause_f lag2=l ; 
DOFOR(k, 19 )  peak_record[ k ]=peak_record( k+1 ] ; 
peak_record[ 19 ]=coordinate[ i  J ; 
if (stability_f lag==l )  stability_check( ) ; 
if (waterfall_f lag==l)  disp_waterf all( ) ; 
) 

else  if (coordinate[ i]<0)  peak_f lag=0; 
) 
)   /*  DOFOR  */ 

BOdel_time=nodel_time+time_int;        /*  INCREMENT  TIME  */ 

/*  CHOOSE  DISPLAY  MODE  AND  EXECUTE  IT  ♦/ 

if (graphics_f lag==l )  display_graphics( ) ; 

if ( ( text_f lag==l ) && ( energy_f lag==l ) )  ( 

printf("\n  Total  Energy  at  time  %lf  is  %lf" ,model_time, energy ) ; 
energy_counter++ ; 
•^    if (energy_counter==20)  ( 

printf ( "\n\nStrike  any  key  to  continue..."); 
while(kbhit( )==0); 
energy_counter=0 ; 
) 
) 

if (phase_f lag!=0)  ( 
phasplot( ) ; 
phase_counter++;      /*  Used  for  Poincare  sections  */ 

) 
}  /*  if(wait_flag  */ 
)  /*  while (stop_f lag. . .  */ 

/*  ***  LEAVE  MAIN  BODY  OF  PROGRAM  HERE  ***  */ 
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/*  ***  EQUATIONS  OF  MOTION.  THE  FOLLOWING  LOOP  CALCULATES  THE 
ACCELERATION  OF  EACH  ELEMENT  AS  THE  FIRST  STEP  OF  THE 
EULER-CROMER  METHOD.  THE  EQUATIONS  CAN  BE  CHANGED  TO 
MODEL  DIFFERENT  PHYSICAL  SYSTEMS  WITH  NO  OTHER  CHANGE 
NEEDED  (ALTHOUGH  IN  ALL  LIKELIHODD  THE  USER_INIT  AND 
FILE  SAVE  ROUTINES  SHOULD  BE  ALTERED  IF  THE  PARAMETERS 
OF  INTEREST  HAVE  CHANGED.   ***  */. 

if (wait_f lag==0)  { 

DOFOR( i , no_pendulums )  (  * 
if(i=0)  { 

acceleration=gamma [ i ] * ( coordinate [ i+1 ] +coordinate [ no_pendulums-l ] 
-2*coordinate[ i ] )-beta*momentum[ i ] 

-(SQR(omegaO[i ] )+2*eta*cos( 2*omega*model_time) )*coordinate[ i ] 
+alpha*CUB(coordinate[i]) ; 
) 
else  if (i==(no_pendulums-l ) )  ( 

acceleration=gamma[ i ]* (coordinate! i-1 }  +  coordinate[ 0  J 
-2*coordinate[ i ] )-beta*moraentum[ i ] 

-(SQR(omegaO[i] )+2*eta*cos(2*omega*model_time) )*coordinate[ i ] 
+alpha*CUB(coordinate[ i ] ) ; 
) 
else  ( 

acceleration=gamma( i ]*(coordinate[ i+1 ]+coordinate[ i-1 ] 
-2*coordinate[ i ] )-beta*momentum( i ] 

-(SQR(omegaO[ i ] )+2*eta*cos(2*omega*model_time) )*coordinate[ i] 
+alpha*CUB(coordinate[ i  J ) ; 
) 
old_momentum[ i ]=momentum[ i ) ; 

momentum!  i ]=momcntum[ 1 ] +acccleration*time_int; 
}   /*  DOFOR  */ 

/*   UPDATE  ACTUAL  ELEMENT  POSITIONS  BASED  ON  THE  NEW  VELOCITIES  */ 

DOFOR (i,no_pendulums)  { 

old_coordinate[ i )=coordinate[ i ] ; 
if (pinned_elements[ i ]==0)  ( 

coordinate!  i  ]=coordinate[  i  )+inonentum[  i  ]*time_int; 
) 

if ((i==chosen_element)ii(file_dump_flag==l) )  (  /*  DUMP  TO  FILE  BUF  */ 
if (counter2++==30)  ( 

time_series_disp[counterl++ J=coordinate( i ] ; 
counter2=0; 
) 

if  (counterl==8000)  stop_f  ile_duinp( ) ; 
)  /*  if((i...  */ 

if (energy_flag==l)  (     /*  CALCULATE  ENERGY  IF  REQUESTED  */ 
energy=energy+0. 5* (SQR( momentum! i ] ) 

+  gamma[ i  J*SQR( (coordinate! i  +  1 J-coordinate[ i  J ) ) 
+SQR( coordinate! i ] ) )-alpha/4*pow(coordinate! i ] , 4 ) ; 

)  /*  i£(encrgy. . . )  */ 

if ((spectrum_f lag==l)66(spectrum_element==i)  /*  UPDATE  SPECTRUM  DATA  */ 
&&( (sample_counter++)==5) )  ( 
if (spectrum_counter<2048)  ( 

spectrum! 2*spectrum_counter ]=coordinatel i ] ; 

spectrum! 2*spectrum_counter+l )=0; 
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_setvideomode(_DEFAULTMODE) ; 

printf  ("PROGRAM  COMPLETE  AT  %lf  "  ,model_time )  ; 
)   /*  MAIN( )  */ 

/****  USER_INIT  ACCEPTS  USER  INPUT  FOR  STARTING  PARAMETERS  AND  INITIAL 
CONDITIONS,  INCLUDING  CHOICE  OF  DATA  FILE  TO  START  FROM.   ALSO, 
VARIOUS  INITIALIZATIONS  ARE  PERFORMED  TO  GET  THE  SYSTEM  GOING  ****/ 

user_init()  { 
FILE  *fq; 
int  c,i,j,k,l,m; 
char  answer[ 1 ] ,f ilename[ 30] ; 

printf ("GENERALIZED  LATTICE  MODEL  PROGRAM  W/  VARIABLE  PARAMETERS\n" ) ; 
printf ("Version  2.1X486  (MS)  \n"); 
printf ("Last  updated  10  JAN  1991\n"); 
printf ("Variant  notes:  Default  1C  is  AM\n"); 
printf ("OmegaO  is  set  equal  to  one  for  all  cases !\n"); 
printf ("Rotating  phase  plane  is  used.\n"); 
printf ("Real  time  FFT  function  is  added. . .\n") ; 
printf ("Energy  monitoring  available  via  *H  in  text  mode"); 
printf ("AL  gives  element  pinning,  NOT  *P. . .\n\n\n" ) ; 

printf ("\nDo  you  want  to  use  a  file  for  initial  conditions  (Y/N)?   "); 
scanf ( "%s" , answer) ; 
Bample_counter=0 ; 

if ((answer[0]==89) | | (answer! 0 )==121 ) )  (   /*  FILE  INPUT  FOR  STARTUP  */ 
printf ("Do  you  want  to  use  old  format  data  file?   "); 
scanf ( "%s" , answer) ; 

printf ("\nEnter  name  of  file  to  be  read:   "); 
scanf ( "%s" , f ilename) ; 
if ((fq=fopen(filename,"r") ) !=NULL)  ( 
mean_gamma=0; 

fscanf ( fq, "%d\n" ,Sno_pendulums) ; 
fscanf (fq,"%lf\n",fcalpha) ; 
fscanf (fq,"% If \n",&beta) ; 
if ( ( answer [0)==89) | | ( answer! 0 )==1 21 ) )  ( 
f scanf ( fq, "%lf \n" , &mean_gamma ) ; 
DOFOR( i ,no_pendulums)  ( 
gamma [ i ] =mean_gamma ; 
omega0[ i ]=1; 
) 
) 

fscanf  (fq,"%lf\n",S.eta)  ; 
fscanf (fq,"%lf\n",iomega) ; 
DOFOR( i ,no_pendulums)  ( 

if ( (answer[0] !=89)&4(answerl0)!=121) )  ( 
fscanf (fq,"*lf   %lf  \lt      %lf\n", 

&omega0[ i ] ,&gamma[ i ] ,fccoordinate[ i ] ,fcmomentum[  i ) ) ; 
ntean_gamma=mean_gamma+gamma  (  i  ] ; 

) 

else  fscanf ( fq, "%lf   *lf \n" ,6coordinate[ i ) , tmomentum[ i ] ) ; 

pinned_elements [ i ) =0 ; 
) 
if ( (answer! 0] ! =B9 )&& ( answer ( 0 ] !=121) ) 

mean_gamma=mean_gamma/no_pendulums; 
f close ( fq) ; 

) 

else  printf ( "Can't  open  file  requested."); 
)  /*  if ( (ans. . .  */ 
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else  {  /*  USER  STARTUP  */ 

printf ( M\nEnter  number  of  pendulums  to  use:   " ) ? 

scanf ( H*d" , &no_pendulums ) ; 

printf ("\nEnter  mode  amplitude:   " ) ; 

ecanf ("%lf ",&mode_amp) ; 

printf ("\nEnter  modulation  amplitude:   "); 

scanf ("%lf " ,&max_amp) ; 

printf ("\nEnter  nonlinear  coefficient  alpha  (+/-  1  ONLY):   "); 

scanf ("%lf",&alpha) ; 

DOFOR (k,no_pendul urns)  ( 

if (alpha<o )  coordinate [ kJ=pow( (-1) ,k)*(mode_amp 
+max_amp*sin(2*PI*k/no_pendulums) ) ; 

else  coordinate[k]«mode_amp+max_amp*sin(2*PI*k/no_pendulums) 

momentum [ k ] =0 ; 

pinned_eleraents [ k J-0 ; 
}   /*  DOFOR  */ 

printf ("\nEnter  coupling  coefficient  gamma:   "); 
scanf ("%lf" ,6mean_gamma) ; 
DOFOR (i ,no_pendulums)  { 

gamma [ i ] =mean_qamma ; 

omega0[i]>l; 
) 

printf ("\nEnter  drive  amplitude  eta:   "); 
scanf ("%lf", Seta); 

printf (n\nEnter  drive  freguency  omega:   "); 
scanf ( "%lf" , iomega ) ; 

printf  (  "\nEnter  dissipation  constant  beta:   •* ); 
scanf  ("%lf",«.  beta); 
)   /*  else  */ 

printf (M\nEnter  time  constant  (multiple  of  period/200):  "); 
scanf (M%lf"#&time_int) ; 
period=2*PI/omega ; 
time_int=time_int*per iod/200 ; 
DOFOR(i,15)  flagsli)=0; 
DOFOR(i,4096)  spectrum[ i ]=0; 
text_flag=l; 
mean_omegaO=l . 0 ; 
spectrum_counter=0 ; 

DINO.02;       /*  CHANGE  THIS  TO  CHANGE  SCALE  OF  DISPLAY  */ 
PINC=2;  /*  CHANGE  THIS  TO  CHANGE  SCALE  OF  PHASE  PLOT  */ 

)  /*  USER_INIT  */ 

/****  GRAPHICS  DISPLAY  ROUTINE  ****/ 

display_graphics( )  ( 
int  c,i,j,k,l,m,n; 
if (no_pendulums<=40)  ( 
f irst_element=0  ; 
DOFOR (k,no_pendulums)  { 
^       n*=0; 

if ((stability_element==k)&&(stability_flag==l))  n=l; 
if ( (pinned_elements[k]==l)&i(disp_color==l ) )  n=2; 
if ( (pinned_elements[k)==l )fc*(disp_color==2) )  n=-l ; 
l-old_coord  i  nate [ k ] /DISPLAY_INCREHENT ; 
_setcolor(0) ; 

_setpixel( (160-5*MIDDLE_ELEMENT  +  5*k) ,( 100+1 )) ; 
l=coordinate[k]/DISPLAY_INCREMENT; 
_setcolor(disp_color+n) ; 
_setpixel((160-5*MIDDLE_ELEMENT  ♦  5*k) ,( 100+1 )) ; 
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)  /*  DOFOR  */ 
)  /*  IF  */ 
else  ( 

l=(no_pendulums-40)/2 ; 
f irst_element=l ; 
DOFOR(k,40)  < 

m=old_coordinate [ 1+k ] /DISPLAY_INCREMENT ; 
n=0; 

if ( (stability_element==(l+k) )&&(stability_f lag==l) )  n=l; 
if  (  (pinned_elements[k]-=l  )fc&(disp_color==l )  )  n=2  ; 
if  ( (pinned_elements[k]==l)&4(disp_color==2) )  n— 1; 
_setcolor(0) ; 

_setpixel((60  +  5*k) , (100+m) ) ; 
in-coordinate  [1+k]  /DISPLAY_INCREMENT; 
_setcolor(disp_color+n) ; 
_setpixel((60  +  5*k) ,  (100+in) )  ; 
)  /*  DOFOR  */ 
)    /*  ELSE  */ 
)  /*  DISPLAY_GRAPHICS  */ 

/****  TEXT  DISPLAY  ROUTINE  ****/ 

display_text( )  ( 

char  message [ 80] ; 

int  c, j ,k,l ,m; 

_setvideomode(_DEFAULTMODE) ; 

printf ("Time  is  :  *lf " ,model_time) ; 

printf("    System  parameters  are:  \n"); 

printf  (•*       Gamma  %lf "  ,mean_gamma )  ; 

print f("      Eta  %lf"  ,eta); 

printf ("       Omega  %lf \n" ,onega ) ; 

printf  ("      Beta  %lf  ••  .beta)  ; 

printf ("      Alpha         %lf\n" .alpha) ; 

printf ("\nThere  are  %d  elements  in  the  system" ,no_pendulums) ; 

printf ( n\n\nPress  any  key  to  continue..."); 

c«getchar( ) ; 

while(kbhit( )==0)  ; 
printf ( "Element      Position      Velocity         Element      Position     Veloci 

DOFOR(j,20)  { 

printf("   %d       %lf  %lf  %d        %lf  % 

j, coord inate[ f irst_element+ j ) , momentum [ f irst_element+ j ] , 
( j+20) .coordinate! first_element+j+ 20 ] • 
momentum[f irst_element+ j  +  20] ) ; 

)  /*  DOFOR  */ 

printf ( "\n\nPress  *G  for  graphics,  'Q  to  quit."); 

c=getchar( ) ; 

while(kbhit()==0); 

while(kbhit()==0); 
}   /*  display_text  */ 

/****  PAUSE  ROUTINE.   THIS  GIVES  USER  CHOICE  OF  SAVING  CURRENT  STATE 

TO  A  DATA  FILE.   ADDITIONALLY,  THE  TIME  STEP  CAN  BE  VARIED  DURING 
A  PAUSE.   ****/ 

pj:ocess_pause  ( )    ( 

int  c,i, j,k,l,mode; 

FILE  *fr; 

char  f ilename[ 30 ] ,  message [80 J; 

_clearscreen(_GCLEARSCREEN) ; 

_setvideomode(_DEFAULTHODE) ; 
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printf("Do  you  wish  to  save  this  state?   "); 
if ((c=getche())==121)  { 

printf ("\nEnter  name  of  file  to  be  written:  "); 
scanf ("%s" , filename) ; 
if ((fr-f open (filename, "w") ) !«NULL)  { 
f printf (fr,"%d\n" ,no_pendulums) ; 
f printf (fr, "%lf\n", alpha ); 
f printf ( f r , "%lf \n" , beta ) ; 
f printf (fr,"*lf\n",eta); 
f printf ( f r , "%lf \n" , omega ) ; 
DOFOR ( i , stabili ty_element ) 

fprintf(fr,"%lf   %lf   %lf  *lf\n", 

omegaO[i ] ,gamma[i] , old_coordinate[ i  ] ,momentum[i] ) ; 
DOFOR(i, (no_pendulums-stability_element) ) 

fprintf (fr,"%lf   %lf   %lf   %lf \n" ,omegaO[ i+stability_element] , 
gamma [i+stability_element ] , coordinate [i+stability_element] , 
momentum [ i+stability_element] ) ; 
fclose(f r) ; 
)  /*  if()  */ 

else  printf ("Failed  to  open  %s\n",f ilename) ; 
)  /*  if  (())  */ 
time_int=time_int*200/period ; 

printf ("\nEnter  new  time  multiple  (old  multiple  is  %lf):   ",time_int); 
scanf  (n%lfn,Sitime_int) ; 
time_int«=time_int*period/200 ; 
if (phase_f lag==0)  _setvideomode(_MRES4C0L0R) ; 
else  _setvideomode(_VRES16C0L0R) ; 
)   /*  process_pause  */ 

/****  START  DUMP  OF  INDIVIDUAL  ELEMENT  TIME  SERIES  TO  DATA  FILE  ****/ 

start_f ile_dump( )  ( 

_clearscreen(_GCLEARSCREEN) ; 

_setvideomode(_DEFAULTMODE) ; 

printf ("Enter  number  of  element  to  be  monitored:  "); 

scanf ("%d" ,&chosen_element) ; 

if (chosen_element>(no_pendulums-l ) ) 

printf ("Out  of  range.   No  file  dump"); 

else  f ile_dump_f lag=l; 

_clearscreen(_GCLEARSCREEN) ; 

if ( phase_f lag==0 ) 

_setvideomode(_MRES4C0L0R) ; 

else 

_setvideomode(_VRES16C0L0R) ; 

counterl=0; 
)  /*  start_f ile_dump  */ 

/****  COMPLETE  DUMP  OF  INDIVIDUAL  ELEMENT  TIME  SERIES  TO  DATA  FILE  ****/ 

stop_f ile_dump( )  ( 

char  filename [30] ; 

int  c; 

FILE  *fs; 

_clearscreen(_GCLEARSCREEN) ; 

_setvideomode ( _DEFAULTMODE ) ; 

printf ("\nEnter  name  of  file  to  be  written:  "); 

scanf ("%s" , filename) ; 

if ((fs»fopen(filename,"w" )) !«NULL)  ( 
■*  DOFOR(c,8000)  f printf ( fs, "tlf\n" , tirae_series_disp(c] ) ; 

_clearscreen(_GCLEARSCREEN) ; 
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counter 1=0; 
if (phase_f lag==0 ) 
_setvideomode(_MRES4COLOR) ; 
else  _setvideomode(_VRES16C0L0R) ; 
)  /*  if(())  */ 
)  /*  Btop_f ile_dump  */ 

/****  CHECK  FOR  SYSTEM  STABILITY.   THIS  IS  A  USER  AID  ONLY,  AND  DOES  NOT 
PERFORM  A  RIGOROUS  STABILITY  CHECK.   IT  SIMPLY  CHECKS  TO  SEE  IF 
ALL  OF  THE  LAST  TWENTY  PEAK  AMPLITUDES  OF  THE  CHOSEN  ELEMENTS  ARE 
WITHIN  ONE  PERCENT  OF  EACH  OTHER  ****/ 

Bonitor_stability( )  ( 

if (stability_flag==0)  { 
stability_f lag=l ; 
disp_color=2; 

_setvideomode(_DEFAULTMODE) ; 

printf( "Enter  number  of  element  to  be  monitored:  "); 
scanf ("%d",fcstability_element); 
stability_restart    ( ) ; 

) 
else  { 

stability_f lag=0; 

disp_color=l; 

) 

if ( phase_f lag==0 ) 
_setvideomode(_MRES4COLOR) ; 
else 

_setvideomode(_VRES16COLOR) ; 
) 

etability_check( )  ( 
int  i,k; 
k=0; 
DOFOR(i,19)  ( 

stable_f lag=l ; 

if ( (peak_record[ i  +  1 ]>( penk_record[ i )*( I « STAB I LITY. JNCRFMntiT ) ) ) | 

(pcnk_record[  Hi  ]<(peak_record(  i  )*(  1  -  STAB1L1TY_1NCREMENT   ))))   ( 
stable_f lag=0 ; 
break; 
)  /*  if  */ 
)    /*  DOFOR  */ 

if (stable_flag==l)  ( 
disp_color=l; 
) 
) 

/****  pin  ANY  INDIVIDUAL  ELEMENT  AT  ZERO  AMPLITUDE.  ****/ 

pin_elements( )  { 
int  ans,i, check; 
_setvideomode(_DEFAULTMODE) ; 
check=0 ; 

printf ( "\nEnter  number  of  element  to  be  pinned:   " ); 
scanf ("%d",ians) ; 

if ( ( ans<0 ) | ( ans> ( no_pendulums-l ) ) )  ( 
check=l; 
) 
if(check==0)  ( 

if (pinned_elements[ans]=»0)  ( 
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pinned_elements[ ans ]=1 ; 

coordinate! ans ]=0; 
momentum [  ans  ]  «=u  ; 

) 

else  pinned_elements[ans]=0; 

)  /*  if  */ 

if (phase_f lag«=0) 
_setvideomode(_MRES4C0LOR) ; 
else 
_setvideomode(_VRES16COLOR) ; 

) 

/****  KICK  ANY  INDIVIDUAL  ELEMENT  ANY  DESIRED  AMOUNT  IN  AMPLITUDE  ****/ 

kick_elements( )  ( 
int  ans, check; 
double  amount; 

_setvideomode(_DEFAULTMODE) ; 
check«0; 

printf (M\nEnter  number  of  element  to  kick:   "); 
scanf ( "%d" , ians ) ; 
if ( (ans<0) | (ans>(no_pendulums-l) ) )  ( 

check=l ; 

) 
if(check==0)  ( 

printf ("\n   Enter  amount  to  kick:   "); 

scanf ("%lf" ,&amount) ; 

coordinate [ ans Incoordinate [ ans ]+amount; 

)  /*  if  */ 
if ( phase_f lag==0 ) 
_setvideomode(_MRES4COLOR) ; 
else 

_setvideomode(_VRES16COLOR) ; 
)  .  ' 

stability_restart( )  { 

int  i  ; 

disp_color=2; 

stable_f lag=0; 

DOFOR(i,20)  peak_record[i]=0; 

peak_record [ 1 5 ] =1 0 ; 
) 

/***#  INITIALIZE  THE  PHASE  PLANE  PLOTTING  SYSTEM.   THIS  ROUTINE  ONTAINS  THE 

USER  CHOICES  FOR  ELEMENTS  TO  MONITOR  AND  DRAWS  THE  AXES  OF  THE  PLOT.  ****/ 

init_phasplot()  ( 
int  c,i,j,k,l; 

char  ans[5] , message [40] ,txt[3] ; 
float  dummy; 

spectrum_f lag=0 ; 

_setvideomode(_DEFAULTMODE) ; 

DOFOR(i,5)  phase_elements[ i ]=0; 

printf ("\nDo  you  want  Poincare  sections?   " ); 

if ( (c=getch( ) )-=121)  phase_f lag=2; 

else  phase_f lag=l ; 

printf ("\nWhich  elements  do  you  wish  to  monitor  (999  to  finish):  "); 

i— 1  ; 

k=0; 


202 


while( (i++!=999)fc&(k<5))  { 
scanf ("%d" ,&i) ; 
phase_elements [k++ ]=i+l ; 
)  /*  while  */ 
DOFOR(i,5)  ( 

if (phase_elements[i ]==1000)  ( 

phase_e lements [ i ]~0 ; 
) 
) 

_setvideomode(_VRES16COLOR) ; 
_clearscreen(_GCLEARSCREEN) ; 
graphics_f lag=0 ; 
_settextcolor(7) ; 
colors [0]=2; 
colors [ 1]=3; 
colors [2 ]=9; 
colors [ 3 ]=14; 
colors [4 ]=13; 
DOFOR(i,120)  ( 

_setpixel(318,4*i); 
) 
DOFOR(i,160)  { 

_setpixel(4*i,238); 
) 
D0F0R(i,6)  ( 

_setpixel ( 319 , 10*SCREEN_CORRECTION_FACTOR+1 
+ ( 80/SCREEN_CORRECTION_FACTOR ) * ( i+ 1 ) ) ; 
) 
D0F0R(i,8)  ( 

_setpixel(80*(i+l) ,239); 
) 

_moveto( 4  80,4  50)  ; 
DOFOR(i,5)  { 

if ( (phase_elements[ i )) !=0)  ( 
c=phase_e lements [ i ]-l ; 
/*       itoa(c,txt,10); 
_outtext(txt);  */ 
printf ("%d   ",c); 
_setcolor ( colors [ i  ]) ; 
D0F0R(j,4)  { 

DOFOR ( k , 8 )  { 

_setpixel(550+10*i+k,10+ j) ; 
) 
) 
) 
)   /*   DOFOR  */ 
printf ("      LAT2X486.C"); 

printf (M\nGamraa:   %lf   Eta:   %lf   Beta:   %lf " ,mean_gamma, eta, beta) ; 
printf ("\nOmega:   %lf   Alpha:   %lf   Time  interval:   %lf" , 

omega , alpha ,time_int ) ; 
dununy=l/(time_int*omega) ; 
phase_int=dummy/l ; 
phase_counter=0 ; 
) 

/****  PHASE  PLANE  PLOTTING  ROUTINE  ****/ 

phasplotO  ( 
int  i,j,k,l; 
double  speed,  position, fast_coordinate,fast_moreentum; 
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if (phaso_f lag==l)  ( 
DOFOU(i,5)  ( 

if  (  (  j^phnr,o_olcmcntn[  i  J-l )  !  — 1 )  ( 

speed=PHASE_INCREMENT*momentum[ j ) ; 
_setcolor(colorsl i ] ) ; 

position=PHASE_INCREMFNT*coordinnto[ j ) ; 
f  ast_coordinate=(  posit  ion*cos( omega*  model  _time) 
-speed*sin(omega*model_time) /omega) 
/( 3*SCREEN_CORRECTION_FACTOR) ; 
f ast_momentum= ( pos  i  t  ion*  s  in ( omega  *mode 1 _t  ime ) 

+speed*cos ( omega*model_t ime ) /omega ) /4 ; 
_setpixel( (320+320*f ast_coordinate) , ( 240+240*f ast_momentum) ) ; 
)  /*  if(())  */ 
)    /*  DOFOR  */ 
)      /*  if  */ 
else  if (phase_f lag==2)  ( 
if  (phase_counter=phase_int)  { 
phase_counter«=0 ; 
DOFOR(i,5)  { 

if  ( ( j=phase_elements[i]-l) !=-l)  ( 

speed=PHASE_INCREMENT*momentum[ j ) ; 
_setcolor( colors [ i  ] ) ; 

position=PHASE_INCREMENT*coordinate[ j ] ; 
fast_coordinate=(position*cos(omega*model_time) 

-speed*sin ( omega*model_t ime ) /omega ) ; 
fast_momentum=(position*sin(omega*model_tine) 

+speed*cos ( omega*model_t ime ) /omega ) /4 ; 
_setpixel( (320+320*fast_coordinate) , ( 240+240*f ast_momentum) ) ; 
)  /*  if(())  */ 
)    /*  DOFOR  */ 
)  /*  if  */ 
)  /*  else  if ( )  */ 
)  /*  phasplot  */ 

/*  FFT  ROUTINE  */ 

compute_spectrum( ) (  /*  Uses  algorithm  from  p.  411  of  Numerical  Recipes  in  C  */ 
int  n,mmax,m, j,istep, i ,nn,temp; 
double  wtemp , wr , wpr , wpi , wi , the ta , tempr , tempi ; 
nn=2048; 
n=nn«l; 

j-i;  .  ' 

DOFOR(i,2048)  ( 

nn=(i&1024)/1024; 

nn=nn+(i&512)/256; 

nn=nn+( i&256 ) /64 ; 

nn=nn+(i&128)/16; 

nn=nn+(i&64)/4; 

nn=nn+(i&32) ; 

nn=nn+(i&16)*4; 

nn=nn+(i&8)*l6; 

nn=nn+(i&4)*64 ; 

nn=nn+(i&2)*256; 

nn=nn+(ifcl)*1024; 

temp2 [ 2*nn )=spectrum[ 2*i ] ; 

temp2[2*nn+l ]=spectrum[2*i+l ] ; 
) 

DOFOR ( i ,4096 )  spectrum[ i ]=temp2 ( i ] ; 
mmax«=2 ; 
while (n>mroax)  ( 
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istep=2*mmax; 
theta=2*PI/mmax; 
wtemp=sin ( 0 . 5*theta ) ; 
wpr=-2 . 0*wtemp*wtemp; 
wpi=sin(theta) ; 
wr«=1.0; 
wi-O.  ef- 
fort m=0;m<(mraax-l )  ;m+=2  )  ( 
for(i=m;i<=n;i+,=istep)  { 
j=i+mmax; 

tempr=wr*spectrum[ j ]-wi*spectrum[ j+1 ] ; 
tempi =wr* spectrum ( j+1 ]+wi*spectrum[ j] ; 
spectrum[ j ]=spectrura( i ]-tempr; 
spectrum[ j+1 ]=spectrum[ i+1 ]-tempi ; 
spectrumf  i ]+=tempr ; 
spectrum [ i+1 ]+=tempi; 
) 

wr= ( wtemp=wr ) *wpr-wi  * wpi+wr ; 
wi=wi*wpr+wtemp*wpi+wi ; 

) 
mmax-=istep; 

) 
) 

/****  SPECTRUM  PLOTTING  ROUTINE  ****/ 

plot_spectrum( )  { 
int  i,j; 

double  freq, magnitude, mag, max; 
int  c; 
FILE  *fp; 
if (dump_spectrum_f lag==l )  { 

printf ( "\nDumping  spectrum  now"); 
f p=fopen( "spectrum. out" , "w" ) ; 
for(i=l;i<2048;i++)  { 

if (time_int!=0)  f req=i/(2048*time_int) ; 
else  printf ("Time_int  was  zero!"); 
magnitude=sqrt(spectrum[ ( 2*i) )*spectrum( ( 2*i ) J 

+spectrum[ ( 2*i+l ) ]*spectrum[ ( 2*i+l ) ] ) ; 
fprintf (fp,n%lf   *lf   \n", freq, magnitude); 

) 

fclose(fp) ; 
dump_spectrum_f lag=0; 

) 

_clearscreen(_GCLEARSCREEN) ; 

max=0 ; 

for(i=l;i<2049;i++)  (  /*  this  loop  calculates  the  max  spectral  component*/ 

magnitude=sqrt( spectrum [ (2*( i ) ) ]*spectrum[ ( 2* ( i ) ) ] 
+spectrum[(2*(i)+l) ]*spectrum[ (2*(i )+l ) j) ; 

if (magnitude>max)  max=magnltude; 
) 

_setcolor(3) ; 

DOFOR(i,640)  _setpixel(i,460) ; 

DOFOR(i,ll)  _setpixel(51*i,461 ) ,_setpixel ( 51 «i ,462) ; 
printf ("Spectrum  for  element  %d      LATTIC2X.C" ,spectrum_element ) ; 
printf ("\nAlpha:  %lf  Beta:  %lf  Gamma:  *lf  Eta:  %lf  Omega:  %lf", 

alpha , beta , mean_gamma , eta , omega ) ; 
printf ("\nMax  amp  =  %lf      ",max); 
freq=512/(2048*time_int) ; 
printf ("Max  freq  shown:  %lf      Time  interval:   tlf " , f req, time_int ) ; 
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_setcolor( 1 ) ; 
DOFOR(i,512)  ( 
j=0; 
magnitude=sqrt( spectrum [ 2*i ]*spectrum( 2*i ] 

+spectrum[2*i+l ]*spectrum[ 2*i+l ] ) ; 
mag=magnitude*400/max; 
while ( j++<mag)  { 

_setpixel(i,460- j) ; 
)  /*  while  */ 
)   /*  DOFOR  */ 
) 

/****  INITIALIZES  WATERFALL  DISPLAY  MODE,  DRAWING  GRID  ****/ 

init_waterfall( )  < 
int  c,i,j,k,l,m; 
char  ans[5]; 

_setvideomode(_VRES16COLOR) ; 
_setcolor(l) ; 

DOFOR( i , 600 )  _setpixel ( 20+i ,  20 ) ; 
DOFOR(i,450)  _setpixel(20,i+20); 
DOFOR(i,61)  ( 

if ( ( i==10 ) | ( i=20 ) | ( i==32 ) | ( i==42 ) | ( i==52 ) )  _setcolor ( 2 ) ; 
DOFOR(j,113)  { 

_setpixel ( 20+10*i , 20+4*  j ) ; 
) 

_setcolor(l) ; 
) 

_setcolor(2) ; 
DOFOR(i,ll)  { 

DOFOR(j,160)  { 

_setpixel( 20+4* j, 16+40* (i  +  1)  ); 
) 
) 

waterf all_counter=0 ; 
color_counter=2 ; 
) 

/****  WATERFALL  DISPLAY  ROUTINE  ****/ 

disp_waterfall( )  { 
int  c,i, j ,k,l; 

if (waterf all_counter++<110)  { 
if (color_counter++"=14)  { 
color_counter=2 ; 
_setcolor(color_counter) ; 
) 
else  ( 

_setcolor(color_counter) ; 
) 

DOFOR (i ,no_pendulums)  ( 
if (waterf all_counter<55)  ( 

_setpixel(20+2*i , 24+8 *water fall  counter-7*wat_inc*( coordinate [ i ] ) ) ; 
) 

else  ( 
_setpixel ( 34 0+2 *i, 24 +8* (waterfall_counter-55) -7 *wat_inc*( coordinate! i ] ) ) » 
) 
) 
) 
else  { 
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